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Chapter  1 

Introduction 


The  project  entitled  “Multiple  Scale  Analysis  of  Damage  and  Texture  Evolution  in  Real  Heteroge¬ 
neous  Materials”  began  in  January  1995.  A  no-cost  extension  was  granted  to  continue  this  project 
tin  December  31,  1998.  During  the  four  year  period  of  this  grant,  substantial  progress  has  been 
made  in  advancing  the  state  of  the  art  in  multiple  scale  modeling  of  damage  in  heterogeneous 
materials.  Research  has  been  conducted  in  a  few  distinct  areas  that  are  delineated  below.  A  list  of 
publications,  acknowledging  this  grant  is  provided  in  chapter  2. 

I.  Voronoi  Cell  Finite  Element  Model  for  Damage  in  Heterogeneous  Microstructures 
(Details  in:  Chapter  3) 

This  work  deals  with  the  evolution  of  damage  in  microstructures  of  reinforced  ductile-matrix  com¬ 
posites,  by  particle  cracking  and  splitting.  A  small  deformation  Voronoi  Cell  finite  element  model 
is  developed,  in  which  each  element  may  consist  of  a  matrix  phase,  an  inclusion  phase  and  a  crack 
phase.  Brittle  inclusions  may  be  of  arbitrary  shapes  and  sizes,  and  may  be  dispersed  non-uniformly 
in  the  matrix.  Damage  initiation  of  inclusions  is  assumed  to  foUow  a  maximum  principal  stress 
theory.  Complete  particle  cracking  or  splitting  is  assumed  at  the  onset  of  damage.  The  model  is 
validated  by  a  few  comparison  studies.  Various  geometric  patterns  are  studied  to  test  the  effective¬ 
ness  of  the  model,  as  well  as  to  understand  the  effect  of  morphology  on  damage  evolution.  Actual 
microstructures  from  optical  micrographs  of  Al-Si-Mg  composite  systems  are  analyzed  and  com¬ 
pared  with  experimentally  observed  results.  Quantitative  characterization  and  statistical  analysis 
is  conducted  to  correlate  morphological  parameters  with  mechanical  response. 

II.  Multi-level  Computational  Model  For  Multiscale  Damage  Analysis  In  Composite 
And  Porous  Materials 

(Details  in:  Chapter  4) 

An  adaptive  multi-level  methodology  is  developed  in  this  work  to  create  a  hierarchy  of  com¬ 
putational  sub-domains  with  varying  resolution  for  multiple  scale  problems.  It  is  intended  to 
concurrently  predict  evolution  of  variables  at  the  structural  and  microstructural  scales,  as  well 
as  to  track  the  incidence  and  propagation  of  microstructural  damage  in  composite  and  porous 
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materials.  The  microstructural  analysis  is  conducted  with  the  Voronoi  cell  finite  element  model 
(VCFEM),  while  a  conventional  displacement  based  FEM  code  executes  the  macroscopic  analysis. 
The  model  introduces  three  levels  in  the  computational  domain  which  include  macro,  macro-micro 
and  microscopic  analysis.  It  differentiates  between  non-critical  and  critical  regions  and  ranges  from 
macroscopic  computations  using  continuum  constitutive  relations  to  zooming  in  at  ‘hotspots’  for 
pure  microscopic  simulations.  Coupling  between  the  scales  in  regions  of  periodic  microstructure 
is  accomplished  through  asymptotic  homogenization.  An  adaptive  process  significantly  increases 
the  efficiency  while  retaining  appropriate  level  of  accuracy  for  each  region.  Numerical  examples 
are  conducted  for  composite  and  porous  materials  with  a  variety  of  microscopic  architectures  to 
demonstrate  the  potential  of  the  model. 

III.  Experimental- Computational  Investigation  Of  Damage  Evolution  In  Discontinu- 
ously  Reinforced  Aluminum  Matrix  Composite 

(Details  in:  Chapter  5) 

This  work  deals  with  a  combined  experimental-computational  approach  to  study  the  evolution 
of  microscopic  damage  to  cause  failure  in  commercial  SiC  particle  reinforced  DRA’s.  Determi¬ 
nation  of  aspects  of  microstructural  geometry  that  are  most  critical  for  damage  nucleation  and 
evolution  forms  a  motivation  for  this  work.  An  interrupted  testing  technique  is  invoked  where 
the  load  is  halted  in  the  material  instability  zone,  following  necking  but  prior  to  fracture.  Sample 
microstructures  in  the  severely  necked  region  are  microscopically  examined  in  3D  using  a  serial 
sectioning  method.  The  micrographs  are  then  stacked  sequentially  on  a  computer  to  reconstruct 
3D  microstructures.  Computer  simulated  equivalent  microstructures  with  elliptical  or  ellipsoidal 
particles  and  cracks  are  constructed  for  enhanced  efficiency,  which  are  followed  by  tessellation  into 
meshes  of  2-D  and  3-D  Voronoi  cells.  Various  characterization  functions  of  geometric  parameters 
are  generated  and  sensitivity  analysis  is  conducted  to  explore  the  influence  of  morphological  pa¬ 
rameters  on  damage.  Micromechanical  modeling  of  2D  micrographs  are  conducted  with  Voronoi 
CeU  Finite  Element  Method  (VCFEM).  Inferrences  on  the  initiation  and  propagation  of  damage  are 
made  from  the  2D  simulations.  Finally,  the  effect  of  size  and  characteristic  lengths  of  representative 
material  element  (RME)  on  the  extent  of  damage  in  the  model  systems  is  investigated. 


IV.  Development  of  the  Voronoi  Cell  Finite  Element  Model  for  Large  Deformation 
Crystal  Plasticity 

A  large  deformation  Voronoi  cell  finite  element  formulation  is  developed  for  modeling  large  elastic- 
plastic  deformation  in  polycrystalline  materials.  The  formulation  is  based  on  as  assumed  stress 
hybrid  finite  element  method  which  makes  independent  assumption  on  stress  rate,  spin  within  the 
domain  of  elements,  and  displacement  or  velocity  assumptions  only  on  the  element  boundary.  The 
finite  element  equations  are  solved  with  a  two  level  Quasi-Newton  iterative  solver.  The  outer  loop 
in  this  iterative  algorithm  is  for  displacements,  and  the  inner  loop  is  for  stress  rates  and  spin, 
with  given  boundary  displacements.  Orders  of  interpolating  polynomials  for  stress  rate  and  spin 
functions  can  be  varied  within  each  element,  independent  of  the  neighbouring  element.  Assuming 
that  slip  is  the  sole  cause  of  plastic  deformation  an  elastic- plastic  constitutive  model  for  crystalline 
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plasticity  is  implemented  into  the  finite  element  algorithm.  A  two-dimensional  crystal  model  with 
three  slip  systems  is  used  as  an  example  model  of  the  single  crystal  that  is  analysed.  The  entire 
crystal  is  modeled  with  one  single  element  in  this  model.  Both  explicit  and  implicit  integration 
procedures  of  the  evolution  equations  are  implemented.  Results  for  elasticity  are  validated  by 
comparison  with  theoretical  solutions  and  results  of  displacement  based  finite  element  codes.  For 
crystal  plasticity,  problems  encountered  modeling  highly  rate  dependent  materials  are  noted  and 
suggestions  to  future  research  to  improve  them  are  made.  A  M.S.  thesis  has  resulted  from  this 
work. 
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Chapter  3 

Voronoi  Cell  FEM  for  Damage  in 
Heterogeneous  Microstructures 


3.1  Introduction 

The  presence  of  fibers  or  particulates  in  composite  microstructures  often  have  adverse  effects  on 
their  failure  properties  like  fracture  toughness,  ductility  and  creep  resistance.  Important  microme¬ 
chanical  phenomena  that  are  responsible  for  deterring  the  overall  properties  include,  fracture  and 
splitting  of  reinforcements,  matrix  failure  and  inclusion-matrix  debonding.  Many  engineering  ma¬ 
terials  exhibit  strong  non-uniformities  in  inter-particle/fiber  spacings,  inclusion  shapes,  volume 
fractions  and  arrangements,  at  the  microstructural  level.  In  addition,  there  are  heterogeneities 
at  larger  length  scales  which  include  local  regions  of  clustering  and  directionality,  often  related 
to  the  fabrication  process.  Failure  characteristics  of  heterogeneous  materials  are  affected  by  mi¬ 
crostructural  mechanisms  that  control  initiation  and  evolution  of  localized  damage  and  cracks. 
These  mechanisms  are  highly  sensitive  to  local  parameters,  such  as  reinforcement  distribution, 
morphology,  size,  interfacial  strength  etc.  Experimental  studies  with  MMCs  [1]  have  established 
that  particles  in  regions  of  clustering  or  preferential  alignment,  have  a  greater  propensity  towards 
fracture,  than  those  in  regions  of  dilute  concentration.  SEM  micrographs  of  damaged  MMCs  show 
that  larger  particles  tend  to  fracture  at  lower  macroscopic  load  levels  due  to  the  existence  of  large 
flaws.  Christman  et.al.  [8]  have  shown  that  local  plastic  flow  is  very  sensitive  to  shape  of  reinforce¬ 
ments.  Lack  of  reliability  of  these  composite  materials  have  inhibited  their  applications  to  high 
performance  load  carrying  engineering  components.  It  is  therefore  important  to  understand  dam¬ 
age  mechanisms  and  fracture  process  for  enhancing  the  level  of  utilization  of  these  material  systems. 

Within  the  framework  of  damage  mechanics,  micromechanical  damage  models  have  been  em¬ 
ployed  to  predict  overall  constitutive  response  by  using  continuum  mechanics  principles  at  the 
microscopic  level  [3,  4].  While  some  of  these  models  [3]  provide  reasonable  predictions  of  overall 
properties  for  a  dilute  distribution  of  damage  entities,  others  [4]  attempt  to  analyze  the  interaction 
effects  between  damage  entities  introduced  by  morphological  characteristics  of  the  microstructure. 
Recently,  novel  approaches  to  integrate  micromechanical  and  computational  approaches  at  the  mi¬ 
croscale  with  phenomenological  approaches  in  the  macroscale  have  also  been  proposed  [5].  While 
many  of  these  methods  can  model  damage  in  brittle  homogeneous  materials,  far  fewer  analytical 
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models  are  available  for  ductile  two  phase  materials.  Small  scale  yielding  solutions  using  asymptotic 
analysis  for  a  single  bi-material  interface  due  to  Shih  &  Asaro  [6]  and  Hutchinson  et.al.  [7]  are 
notable  exceptions. 

Evolving  damage  in  heterogeneous  media  with  a  mixture  of  ductile  and  brittle  constituents  have 
been  numerically  modeled  using  Unit  Cell  methods.  These  methods  assume  that  the  material  is 
constituted  of  periodic  repetition  of  unit  cells,  identified  as  representative  volume  elements  (RVE) 
of  the  microstructure.  Displacement  based  finite  element  analysis  is  used  to  analyze  the  RVE  in 
order  to  predict  the  onset  and  growth  of  evolving  damage.  Notable  among  these  are  the  finite 
element  simulations  by  Needleman  [11,  9],  Tvergaard  [12],  Bao  [13],  Horn  [14]  ,  Sugimura  et.al. 
[13]  and  Finot  et.al.  [15].  In  [12,  13,  14,  13],  simple  microstructures  with  pre-existing  damaged 
heterogeneities  are  considered.  In  [15]  a  finite  element  mesh  which  allows  for  crack  growth  by  ele¬ 
ment  separation  is  used  to  simulate  microscale  particle  and  matrix  cracking.  While  these  models 
provide  valuable  insights  into  the  microstructural  damage  processes,  the  simple  morphologies  ide¬ 
alize  actual  microstructures  for  many  engineering  materials.  Consideration  of  simple  RVE’s  bear 
little  relationship  to  the  actual  stereographic  features,  and  have  limited  them  to  the  assumption 
that  all  particles  or  particle/matrix  interfaces  are  damaged  simultaneously.  To  circumvent  these 
deficiencies,  Suresh  and  coworkers  [16,  8],  McHugh  et.  al.  [17]  among  others,  have  made  novel 
progresses  in  computational  modeling  of  discontinuously  reinforced  materials  with  random  spatial 
dispersion.  However,  a  very  high  resolution  of  finite  element  mesh  is  required  even  for  undamaged 
heterogeneous  media,  and  enormous  computational  efforts  are  required  to  capture  failure  by  these 
models. 

The  microstructure  based  Voronoi  Cell  Finite  Element  Model  (VCFEM)  developed  by  Ghosh 
et.al.  [34,  35],  has  shown  a  significant  promise  in  this  regard.  It  can  overcome  the  large  com¬ 
putational  requirements  of  conventional  finite  element  methods,  by  combining  concepts  of  hybrid 
finite  elements  with  characteristics  of  micromechanics.  The  VCFEM  mesh  naturally  evolves  from 
the  microstructure  by  Dirichlet  tessellation  to  generate  a  network  of  multi- sided  Voronoi  polygons. 
Each  Voronoi  cell  represents  a  basic  structural  element  containing  one  second  phase  inclusion  at 
most  (see  [19]  for  details),  and  the  analysis  needs  no  further  discretization  leading  to  drastically 
reduced  efforts  in  generating  the  microstructural  mesh.  Additionally,  computational  efficiency  is 
greatly  enhanced  due  to  Voronoi  cell  elements  being  considerably  larger  than  conventional  unit  ceU 
finite  elements,  with  reduced  degrees  of  freedom. 

The  evolution  of  damage  by  particle  cracking  or  splitting,  in  particle  reinforced  ductile  matrix 
microstructures,  is  analyzed  in  this  work  by  a  Voronoi  ceU  finite  element  model.  No  matrix  cracking 
is  allowed  in  this  work.  Each  Voronoi  ceU  element  may  consist  of  a  matrix  phase,  an  inclusion  phase 
and  a  crack  phase.  The  inclusions  are  brittle,  of  arbitrary  shapes  and  sizes,  and  may  be  dispersed 
non-uniformly  in  the  matrix.  Damage  initiation  is  assumed  to  foUow  a  maximum  principal  stress 
theory  or  Rankine  criterion.  Complete  particle  cracking  or  splitting  is  assumed  at  the  onset  of 
damage.  Different  geometric  patterns  are  studied  to  test  the  effectiveness  of  the  model,  as  weU 
as  to  understand  the  effect  of  morphology  on  damage  evolution.  Actual  microstructures  from 
optical  micrographs  of  Al-Si-Mg  composite  systems  are  analyzed  and  compared  with  experimentally 
observed  results.  Quantitative  charact  erization  and  statistical  analysis  is  conducted  to  correlate 
morphological  parameters  with  mechanical  response. 
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3.2  Voronoi  Cell  FEM  with  Particle  Fracture 


The  Voronoi  cell  finite  element  model  has  been  developed  for  undamaged  composite  and  porous 
materials  in  [34,  35]  using  the  assumed  stress  hybrid  formulation.  The  formulation  is  extended 
to  accommodate  damage  evolution  in  the  form  of  particle  cracking  or  splitting.  It  is  assumed 
that  in  particle  cracking,  the  crack  is  completely  contained  within  the  inclusion,  while  for  particle 
splitting  it’s  tip  extends  nominally  into  the  matrix.  The  crack  in  a  fractured  particle  is  realized  as 
an  elliptical  void  with  a  high  aspect  ratio  (~10  -100),  implying  a  blunt  crack.  Each  Voronoi  ceU 
element  is  amenable  to  change  in  topology  from  two  constituent  phases  (matrix  and  inclusions)  in 
undamaged  cells,  to  three  phases  (matrix,  inclusion  and  crack)  in  damaged  cells.  Complete  particle 
cracking  or  splitting  is  assumed  to  occur  at  the  very  onset  of  damage,  and  thus  the  problem  of  crack 
propagation  within  each  inclusion  is  avoided.  This  assumption  is  justifiable  from  the  consideration 
that  for  the  multitude  of  inclusions  analyzed,  crack  propagation  in  each  inclusion  would  make  the 
problem  inordinately  large.  Additionally,  experimental  observations  indicate  rapid  transition  from 
crack  initiation  to  complete  cracking/splitting. 


3.2.1  Voronoi  cell  element  formulation  for  damage 

Consider  a  typical  representative  material  element  (RME)  consisting  of  N  undamaged  and/or 
damaged  particles,  that  are  contained  in  each  of  the  N  Voronoi  cell  elements  as  shown  in  figure 
3.5(a).  The  assumed  stress  hybrid  formulation  in  the  Voronoi  ceU  finite  element  method  (VCFEM) 
requires  independent  assumptions  of  an  equilibriated  stress  field  (c)  in  the  interior  of  each  element 
fte,  and  compatible  displacement  fields  u  on  the  element  boundary  dfte,  on  the  matrix-inclusion 
interface  dClc  and  u''  on  the  crack  boundary  d^cr-  la  an  incremental  formulation  for  elasto- 
plasticity,  the  incremental  two  field  (<r  —  u)  hybrid  variational  formulation  introduces  an  element 
energy  functional. 


n^( A<r,  Au)  =  —  /  A5(«r,  A<r)  dfl  —  f  e:A<rdQ  + 

fie 

f  {<T  +  Acr)  •  n®  •  (u  -|-  Au)  dO,—  (  (t  -|-  At)  •  (u  -t-  Au)  dT 
JdQe  JVtm 

-  f  (o-”^  +  Ao-’"  -  -  Atr^)  •  •  (u'  +  Au')  dCl 

JdQc 

-  /  +  Ao-")  ■  .  (u"  -h  Au")  dU 

JdQcr 


(3.1) 


where  AB  is  the  increment  of  complimentary  energy  density.  Variables  (<t,u)  correspond  to  values 
at  the  beginning  of  an  increment,  while  variables  (A<r,  Au)  are  the  corresponding  increments  in  a 
load  increment  or  step.  Outward  normals  on  dilg,  dSlc  and  dQcr  are  denoted  by  n®,  n®  and  n®’’ 
respectively.  Superscripts  m,  c  and  cr  are  associated  with  the  matrix,  inclusion  and  crack  phases 
respectively  in  each  Voronoi  ceU  element.  The  total  energy  for  the  entire  RME  of  N  Voronoi  ceUs  is 
obtained  as  11^  =  •  Setting  the  first  variation  of  11^  in  equation  (3.1)  with  respect  to  stress 

increments  A<r  to  zero  yields  the  element  compatibility  as  the  Euler  equation,  while  setting  the 
first  variations  of  11^  with  respect  to  the  independent  boundary  displacements  Au,  Au'  and  Au" 
to  zero,  yield  the  inter-element  traction  reciprocity  or  element  boundar  traction,  interface  traction 
reciprocity  and  zero  traction  on  crack  boundary  respectively.  Equilibriated  stress  increments  Acr, 
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compatible  displacement  fields  Au,  Au'  and  Au",  the  stress-strain  relationships  (§^§^  =  Ae), 
along  with  the  Euler  equations  completely  define  the  incremental  problem  for  a  heterogeneous 


RME. 


3.2.2  Element  assumptions 

Independent  assumptions  on  stress  increments  Acr  are  made  in  the  matrix  and  inclusion  phases 
in  each  element,  thus  allowing  stress  discontinuities  across  the  interface.  In  two-dimensional  anal¬ 
ysis,  the  Airy’s  stress  function  ^(x,y)  is  usually  convenient  in  deriving  equilibriated  stress  fields. 
Components  of  Ac  are  expressed  in  terms  of  $  as: 

A  _  _  A  _  A  _  _ 

""  -  0y2  ’  -  Q^Qy 

Incorporation  of  key  features  of  micromechanics  in  the  choice  of  stress  functions  significantly  en¬ 
hances  computational  efficiency.  Moorthy  and  Ghosh  [34,  19]  have  introduced  a  decomposition  of 
the  matrix  and  inclusion  stress  functions  into  (a)  purely  polynomial  functions  ^poly  (^) 

reciprocal  functions  ^rec  ^reo  for  elements  with  matrix,  inclusion  and  crack  phases. 

Mathematically,  the  stress  functions  for  the  matrix  and  inclusion  phases  are  constructed  as: 


■ar  —  ^poly  ~  ^rec  ~  ^rec 

-f 

(3.3) 

In  the  above  equation,  the  purely  polynomial  part  of  the  stress  functions  account  for  the  far 

field  tractions  on  the  element  boundary  dClg  and  on  the  interface  are  expressed  as: 

(3-4) 

where  (^,  rj)  are  the  scaled  local  coordinates  with  origin  at  the  element  centroid  (xc,  yc),  and  may  be 
written  as  ^  =  {x—Xc)/L,  t)  =  {y—yc)lL.  i  is  a  scaling  parameter  (=  y/max{x  —  x^)  max{y  —  yc)  V(x,y)  €  d^. 
The  use  of  the  local  coordinates  (£,  77)  prevents  numerical  inaccuracies  in  due  to  high  expo¬ 
nents  of  (x,  y),  and  thus  avoid  ill-conditioning  of  the  element  stiffness  matrices.  The  reciprocal 
terms  facilitate  stress  concentration  near  the  interface  and  crack  boundary, 

accounting  for  the  shape  of  the  inclusion  and  crack.  They  also  help  satisfy  traction  reciprocity 
(zero  traction  for  the  crack)  at  the  interfaces  and  dO.cr,  as  well  as  decay  at  large  distances 
from  these  interfaces.  The  matrix  reciprocal  function  is  constructed  from  a  transformed  radial 
coordinate  /,  that  is  generated  by  either  a  Schwarz- Christoffel  conformal  transformation  (for  ellip¬ 
tical  heterogeneities)  [20],  or  by  a  Fourier  series  transformation  of  the  interface  dQ.c  (for  arbitrary 
shapes)  [34].  The  radial  distance  /  satisfies  the  conditions  /  ^  00  as  (x,y)  — »•  00,  and  /  =  1  on 
the  interface  d^c-  In  the  expression  for  shape  effects  are  dominant  near  d0.c  and  vanish  in 
the  far-field. 

^rec  =  (3-5) 

P>9  »  •' 

At  the  interface  (/  =  1),  coefficients  A/3^,-  in  equation  (3.5)  impart  flexibility  to  the  polynomial  co¬ 
efficients  A/?^  for  matching  traction  conditions.  Finally,  the  terms  are  contributions 
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to  the  matrix  and  inclusion  stress  functions  due  to  the  crack. 

EfVD7F?fe) 

p,g  i  Jcr 

(3.6) 

P,9  t 

The  inclusion  crack  is  assumed  to  be  of  elliptical  shape  with  a  high  aspect  ratio.  Consequently, 
the  crack  boundary  d^la-  is  parametrically  represented  through  a  conformal  mapping  of  the  ellipse 
fcrix,y)  =  1.  fcr  represents  a  parameterized  radial  coordinate  with  the  property  /cr  — >  oo  as 
i^^y)  00.  The  reciprocal  terms  in  facilitate  zero  traction  condition  on  the  crack 

boundary  ddcr-  The  same  terms  in  provide  asymptotic  stress  gradients  near  the  crack  tip  in 

the  matrix.  Stress  increments  may  be  derived  by  substituting  ^  functions  in  the  equation  (3.2)  in 
the  form  {Aff”*}  =  [P”^]{A;3'"}  for  the  matrix  and  {A<r^}  =  [P^]{A/3‘^}  for  the  inclusion.  AH  of 
the  stress  coefficients  and  are  a-priori  unknown  and  are  solved  by  setting  the  first 

variation  of  the  element  energy  functional  (3.1)  with  respect  to  the  stresses  to  zero.  Compatible 
displacement  increments  are  generated  on  each  of  the  boundaries/interfaces  dCle,  dQc  and  dQcr  by 
interpolation  in  terms  of  generalized  nodal  displacements  as  , 

{Au}  =  [L']{Aq}  ,  {Au'}  =  [L^]{Aq'}  ,  {Au"}  =  [L-]{Aq"}  (3.7) 

where  {Aq},  {Aq'}  and  {Aq"}  are  the  nodal  displacement  increment  vectors,  and  [L®],  [L*^]  and 
are  the  corresponding  interpolation  matrices.  In  general,  linear  forms  of  [L]  are  computation¬ 
ally  efficient.  However  for  the  crack  boundary,  discontinuous  normals  at  the  nodes  may  degrade 
the  solution  and  hence  a  quadratic  interpolation  is  implemented.  Details  of  the  solution  process 
are  provided  in  [34,  19,  21]. 

3.2.3  Constitutive  relations  and  particle  cracking  criterion 

The  reinforcing  phase  of  particles  are  assumed  to  be  brittle  and  are  modeled  as  linear  elastic 
materials.  The  matrix  material  on  the  hand  is  assumed  to  be  ductile,  and  is  modeled  by  small 
deformation  elasto-plasticity  relations  using  associated  J2  flow  theory  with  isotropic  hardening.  For 
the  brittle  particulate  materials,  microstructural  damage  initiation  is  assumed  to  be  governed  by 
a  maximum  principal  stress  based  criterion,  also  known  as  the  Rankine  criterion.  In  this  criterion, 
a  crack  is  initiated  when  the  maximum  principal  stress  in  tension  exceeds  a  critical  fracture  stress 
(Tct  at  a  point.  In  the  computational  procedure,  complete  particle  cracking  or  splitting  is  assumed 
to  occur  in  the  form  of  an  elliptical  void,  as  soon  as  the  principal  tensile  stress  reaches  acr-  In 
the  case  of  particle  cracking,  the  crack  tip  coincides  with  the  interface  and  is  completely  contained 
in  the  particle,  while  it  extends  nominally  into  the  matrix  for  particle  splitting.  A  parameter 
^crack  =  Inclusion  Dmension  distinguishes  between  Complete  cracking  and  splitting  of  inclusions.  A 
fuUy  cracked  inclusion  corresponds  to  a  value  dcrack  =  Ij  for  which  the  crack  terminates  at  the 
inclusion-matrix  interface,  whereas  splitting  is  represented  by  dcrack  =  1.004  for  which  the  crack 
tip  has  moved  slightly  into  the  matrix.  In  the  incremental  computational  procedure,  more  than 
one  point  may  exceed  the  critical  acr  value  during  increment.  The  location  of  a  single  crack  is 


_ 

^rec  ““ 
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determined  by  a  weighted  averaging  method  as: 


^damage  = 


ydamage  — 


Ey 


cr 


Egf(^.y) 
fcr 


V  [aKa:,  y)  >  (Ter] 


(3.8) 


where  (Tj(a:,  j/)  corresponds  to  aU  values  of  maximum  tensile  principal  stress  larger  than  <7^  in  the 
particle.  The  crack  is  oriented  at  right  angles  to  the  principal  stress  directions  at  {x damages  y damage) 
and  extends  to  the  interface  on  both  sides. 


Critical  Fracture  Stress 


Various  experimental  studies  on  metal  matrix  composites  [6,  22,  24],  suggest  that  the  critical  stress 
O’er  for  particle  fracturing  is  not  only  material  dependent,  but  is  also  influenced  by  the  particle 
size  due  to  the  existence  of  microcracks.  Micrographs  of  damaged  composites  indicate  that  larger 
particles  tend  to  fracture  at  lower  load  levels  than  smaller  particles.  To  account  for  the  size  eflFect  in 
particle  cracking,  and  hence  flaw  size  and  distribution,  two  alternative  approaches  are  considered. 
These  criteria  have  been  discussed  in  Curtin  [25],  Kiser  et.al.  [24].  The  first  is  a  fracture  mechanics 
based  criterion,  in  which  particles  are  assumed  to  contain  flaws  and  the  critical  stress  to  fracture 
is  determined  from  mode-I  fast  fracture  of  these  flaws.  In  this  criterion,  an  initial  particle  flaw 
size  c  is  assumed  to  be  a  fraction  of  a  charcicteristic  length  D,  and  is  expressed  as.  c  =  eD.  The 
characteristic  length  is  considered  to  be  the  diameter  of  an  equivalent  circle  or  D  =  where 

A  is  the  particle  area.  The  factor  e  is  determined  from  experimental  observations,  and  a  value 
~  5  —  15%  is  found  to  be  suitable  in  this  study.  For  mode-I  fracture,  the  critical  load  to  fracture 
(Tct  is  thus  related  to  the  fracture  toughness  Kjc  through  the  relation: 


(Ter 


Kic  ^  Kjc 
y/^  Vtt  eD 


(3.9) 


Larger  particles  with  large  initial  flaws  will  fracture  at  smaller  critical  stresses  by  this  relation.  The 
second  criterion  uses  statistical  functions  to  correlate  particle  size,  stress  levels  and  failure.  It  is 
based  on  a  WeibuU  distribution,  in  which  the  probability  of  particle  fracture  P/(A,cr)  is  related  to 
the  particle  volume  (area  in  2-D)  A  and  the  maximum  principal  stress  aj  as: 


P/(A,  <7)=1  — e  ^(<Tcr)  (3.10) 

where  (7^-  and  m  are  two  material  parameters  in  the  WeibuU  distribution.  The  probabiUty  of  damage 
in  this  model,  increases  with  larger  particles  at  larger  stress  levels.  The  WeibuU  parameters  Ocr  and 
m  may  be  calculated  by  correlating  geometric  features  and  simulated  stresses  with  experimental 
observations,  as  discussed  in  the  section  on  numerical  examples. 


3.3  Validation  Examples  for  Voronoi  Cell  FEM 

The  accuracy  and  efficiency  of  the  Voronoi  ceU  finite  element  model  in  stress  analysis  of  hetero¬ 
geneous  materials  with  particle  cracking  has  been  extensively  verified  by  comparison  with  results 
of  analyses  conducted  with  conventional  finite  element  packages  as  weU  as  with  pubUshed  results 
in  the  Uterature.  Several  comparison  studies  have  been  made  with  this  model,  some  of  which  are 
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described  in  [21, 19].  Only  one  such  example  studies  with  simple,  uniform  distribution  is  described 
here.  To  establish  an  aspect  ratio  for  the  elliptical  cracks,  a  numerical  experiment  was  conducted 
with  various  aspect  ratios  viz.  (|  =  3,5,10,25,20  and  100).  From  this  study,  a  ratio  f  =  10  was 
found  to  be  most  desirable  from  an  efficiency  and  accuracy  point  of  view.  While  both  plane  stress 
and  plane  strain  problems  have  been  solved,  the  numerical  examples  presented  in  this  work  are 
only  for  the  plane  strain  assumption. 

3.3.1  Square  diagonal  packing  with  existing  crack 

This  comparison  problem,  studied  by  Finot  et.al.  [15]  with  finite  deformation  kinematics,  involves 
stress  analysis  of  a  square-diagonally  packed  composite  microstructure  with  pre-cracked  inclusions. 
The  representative  material  element  (RME)  consists  of  two  square  SiC  inclusions  (volume  fraction 
Vj  =  20%)  in  an  Al-3.5%Cu  matrix  as  shown  in  figure  (3.3.1).  The  elastic  properties  of  SiC  par¬ 
ticles  are  assumed  to  be:  Young’s  Modulus  E  =  450  GPa  and  Poisson’s  Ratio  u  =  0.2.  The 
elastic-plastic  properties  for  the  Al-3.5%  Cu  alloy  matrix  are  taken  as:  Young’s  Modulus  E  =  72 
GPa,  Poisson’s  Ratio  v  =  0.32;  Post  yield  behavior  (Power  law  hardening)  Om  =  cro(€^/eo  -1- 1)^, 
with  (To  =  175  MPa,  and  N  =  0.2.  Different  degrees  of  pre-existing  damage,  e.g.  0%  damage  with 
two  intact  inclusions,  50%  damage  with  cracked  inclusion,  and  100%  damage  with  both  inclusions 
cracked,  are  assumed  in  accordance  with  those  used  in  [15].  For  0%  damage,  the  matrix  stress 
function  in  equation  (3.3)  consists  of  61  terms,  with  25  polynomial  terms  (p  +  9  =  2.. 6 
in  equation  3.4)  and  36  reciprocal  terms  (*  =  1..3,p-l-  q  =  2. .6  in  equation  3.5).  The  cor¬ 
responding  inclusion  stress  function  in  equation  (3.3)  consists  of  25  polynomial  terms 
{p  +  q  =  2. .6  in  equation  3.5).  The  function  /  in  equation  (3.5)  for  is  created  by  a  Fourier 
series  transformation  of  the  square  interface,  as  described  in  [34].  For  damaged  microstructures, 
36  additional  terms  in  the  form  of  and  are  appended  to  the  stress  functions  and 
(i  =  1..3,p-f  q  =  2. .4  in  equation  3.6). 

The  VCFEM  simulation  of  the  RME  is  executed  for  upto  a  vertical  applied  strain  of  2%  as 
shown  in  figure  (3.3.1).  In  figure  (3.3.1)  the  macroscopic  stress-strain  responses,  calculated  by 
taking  volumetric  averages  of  microscopic  variables,  are  compared  with  results  in  [15],  and  excellent 
agreement  is  recorded.  The  stress  capacity  of  the  RME  reduces  considerably  with  transition  from 
particle  cracking  to  particle  splitting.  As  the  crack  propagates  into  the  matrix  due  to  splitting, 
the  damaged  inclusions  cease  to  carry  significant  load.  The  major  load  now  shifts  to  the  matrix 
material  and  the  remaining  undamaged  inclusions.  Contour  plots  of  the  effective  plastic  strains 
for  the  cracked  and  split  microstructures  at  Cyy  =  2%  are  presented  in  figures  (3.3.1).  The  matrix 
regions  vertically  adjacent  to  the  split  inclusions  have  considerably  less  plastic  strains  than  those 
adjacent  to  cracked  inclusions,  due  to  much  lower  stresses  caused  by  splitting.  Also  in  the  case 
of  splitting,  a  considerably  larger  plastic  strain  accumulates  near  the  crack  tip.  The  plastic  strain 
flows  in  the  form  of  ligaments  from  one  crack  tip  to  the  next,  causing  bands  of  strain  localization. 
Similar  observations  have  also  been  made  in  [15]  for  axisymmetric  inclusions. 


3.4  Damage  in  Non-uniform  MicrostructureS 

Examples  in  the  previous  section  consider  pre-existing  damage,  and  thus  do  not  involve  crack 
initiation  and  change  of  element  topology.  The  present  section  deals  with  continuously  evolving 
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Figure  3.1:  Representation  and  VCFEM  mesh  for  RME  with  V/  =  20%  cracked  square  inclusions 
(a)  0%  jg^g.ge  (b)  50%  damage  ^c)  100%  damage. 
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Figure 


Macroscopic  Tensile  Strain  (%) 

3.2:  Macroscopic  stress-strain  response  for  RMEs  with  Vf  =  20%  cracked  square  inclusions. 
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microstructural  topology  through  the  onset  and  evolution  of  particle  cracking  in  more  complex, 
computer  simulated  and  real  microstructures.  For  undamaged  Voronoi  elements,  the  matrix  stress 
function  in  equation  (3.3)  consists  of  34  terms  with  25  polynomial  terms  (P  +  9  =  2. .6 
in  equation  3.4)  and  9  reciprocal  terms  (z  =  1..3,p  =  0..2,  q  =  2  —  p  in  equation  3.5).  The 
inclusion  stress  function  consists  of  25  polynomial  terms  (p  +  q  =  2. .6  in  equation  3.5). 

Fewer  reciprocal  terms  are  used  than  in  the  previous  example  due  to  the  smooth  interface  for 
circular  inclusions.  For  cracked  elements  though,  the  additional  reciprocal  terms  are  the  same  as 
those  in  the  previous  example,  i.e.  {i  =  1..3,p  +  q  =  2. .4  in  equation  3.6). 

3.4.1  Computer  simulated  microstructures 

The  effect  of  particle  clustering  on  damage  evolution  is  studied  with  two  computer  generated 
microstructural  distributions  as  follows. 

(a)  A  hard  core  distribution:  which  is  generated  as  a  variant  of  a  pure  random  Poisson  pattern 
through  the  imposition  of  two  constraints,  namely  (a)  no  two  inclusions  are  allowed  to  overlap,  and 

(b)  all  inclusions  are  completely  contained  within  the  region  window. 

(b)  A  single  cluster  hard  core  distribution  . 

(c)  Triple  cluster  hard  core  model,  where  clusters  are  characterized  by  reduced  average  inclusion 
separation  distance  in  an  otherwise  hard-core  region. 

Each  RME  consists  of  50  equi-sized  circular  Si  particles  dispersed  in  an  Al-Si-Mg  alloy  matrix,  and 
constituting  a  20%  volume  fraction.  Pertinent  dimensions  are  :  RME  size=200/x  x  200/i,  particle 
diameter  =14.2/i,  cluster  diameter  in  (b)  =  33.38//  and  cluster  diameter  in  (c)=  25.78//.  Details  of 
the  generation  process  are  described  in  [27]. 

Characterization 

Statistical  functions  of  geometric  descriptors,  which  can  discriminate  between  patterns,  as  discussed 
in  [27,  28],  are  considered.  Figures  3.4.1(a  and  b)  show  the  cumulative  distribution  function  F{A) 
and  the  probability  density  function  f{A)  of  the  local  area  fraction,  measured  as  the  ratio  of  the 
particle  size  to  the  area  of  the  associated  Voronoi  cell.  The  range  of  A  for  hard  core  pattern  is 
significantly  shorter  than  for  clustered  patterns  and  thus  the  difference  in  F(A)  increases  with 
increasing  area  fraction.  The  high  spike  in  f(A)  for  the  hard  core  pattern  is  a  consequence  of  the 
steep  gradients  due  to  pronounced  uniformity  in  local  area  fraction,  and  the  intensity  of  these  spikes 
diminishes  with  clustering.  The  cumulative  distribution  function  F{d)  and  density  distribution 
functions  f(d)  for  center  to  center  nearest  neighbor  distances  are  plotted  in  figures  3.4.1(c  and  d). 
The  longer  plateaus  in  F(d)  and  the  corresponding  zeros  in  f{d)  for  clustered  patterns  are  for  the 
distances  for  which  a  near  neighbor  does  not  exist.  Spikes  in  f(d)  are  much  more  pronounced  for 
the  hard-core  distribution  due  to  large  number  of  neighbors  at  nearly  similar  distances.  The  lowest 
d  values  are  much  smaller  for  the  clustered  patterns.  Second  order  intensity  function  K{r),  defined 
as  the  number  of  additional  points  expected  to  lie  within  a  distance  r  of  an  arbitrarily  located 
point  divided  by  overall  the  point  density,  is  an  informative  descriptor  and  has  been  discussed 
in  [22,  27,  28]).  Additionally,  the  pair  distribution  function  q(r)  =  corresponds  to  the 

probability  g{r)dr  of  an  finding  an  additional  point  within  a  circle  of  radius  dr  and  centered  at 
r.  The  two  functions  are  plotted  for  the  patterns  in  figure  (3.4. le  and  f)  and  compared  with 
a  pure  Poisson  process  for  which  K(r)  =  zrr^  and  g(r)  =  1.  With  increased  in  clustering,  K(r) 
deviates  from  that  for  the  Poisson  process.  The  peaks  in  q(r)  are  more  pronounced  for  the  hardcore 
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distribution  indicating  a  greater  likelihood  of  encountering  additional  particles  at  lower  radii  for 
this  volume  fraction. 

Damage  Simulation 

Both  complete  particle  cracking  and  particle  splitting  are  analyzed  in  this  example,  with  initially 
undamaged  particles.  Constituent  material  properties  are  as  follows; 

For  the  Al-Si-Mg  matrix:  Young’s  modulus  E  =  69  GPa,  Poisson’s  ratio  u  =  0.33,  and  the  post 
yield  elastic-plastic  behavior  is  obtained  from  data  on  B-treatment  matrix  material  in  figure  (5-1) 
(page  137)  of  Hunt  [22].  For  the  Si  inclusions:  Young’s  modulus  E  =  161  GPa,  Poisson’s  ratio 
u  =  0.2.  AH  particles  are  of  identical  shape  and  size,  and  therefore  a  constant  critical  stress  to 
failure  <7^  =  300  MPa  is  assumed.  The  RME’s  are  subjected  to  a  macroscopic  horizontal  tensile 
strain  that  increases  from  0  to  a  maximum  of  l^x  =  2%.  Evolving  damaged  configurations  by 
particle  splitting  are  shown  in  figures  (3.5a-f).  For  the  hard  core  pattern,  the  first  set  of  particles 
crack  at  Ixx  =  0.6%.  Particle  cracking  occurs  at  random  locations  with  increased  loading  upto 
a  strain  of  =  1.6%,  after  which  no  additional  cracking  is  noticed.  Large  plastic  strains  occur 
in  regions  near  the  crack  tips.  However  due  to  the  lack  of  close  proximity  of  cracked  particles, 
plastic  strain  regions  were  not  observed  to  propagate  in  any  preferred  direction.  In  contrast,  for 
the  triple  cluster  microstructure  in  figure  3.5d-f,  the  first  set  of  particles  crack  and  spht  within  the 
cluster  near  the  top-right  corner  at  e^x  =  0.4%.  With  increased  loading,  high  stress  concentration 
at  the  crack  tips  lead  to  progressive  cracking  of  other  particles  inside  this  cluster  before  particles  in 
other  clusters  begin  to  crack.  However  if  only  particle  cracking  is  allowed,  particles  in  both  the  top 
right  and  bottom  center  clusters  begin  to  crack  almost  concurrently.  Particles  in  the  third  cluster 
remain  rel  atively  undamaged  during  the  entire  process.  At  the  final  strain  Ixx  =  2.0%,  most 
particles  in  the  two  clusters  are  split  while  damage  in  the  third  cluster  has  just  begun.  During  the 
initial  stages  of  deformation  and  splitting,  localized  plastic  straining  occurs  in  the  top  right  cluster 
which  propagates  from  one  crack  tip  to  the  next  within  the  cluster  by  linking.  Plastic  straining  in 
other  clusters  are  less  pronounced  during  this  period.  With  subsequent  particles  splitting,  plastic 
straining  intensifies  in  bottom-center  cluster  and  eventually  links  up  with  strained  regions  in  the 
top  cluster.  High  strain  regions  are  much  more  diffused  in  the  case  of  particle  cracking  and  occur 
at  higher  macroscopic  stresses  compared  to  the  particle  splitting  case. 

Macroscopic  stress-strain  responses  of  the  hard  core  and  triple  clustered  RMEs  are  illustrated 
in  figure  (3.4.1).  Abrupt  drops  due  to  particle  cracking  are  smoothed  in  this  figure.  The  stress  level 
for  the  hard  core  pattern  with  particle  splitting  continues  to  drop  throughout  the  loading.  For  the 
clustered  microstructure,  drops  are  higher  in  the  initial  stages  due  to  rapid  failure  in  the  clusters, 
followed  by  increase  in  the  stress  levels  due  to  matrix  hardening.  In  general  the  RME  with  cracked 
particles  projects  a  considerably  stiffer  behavior  when  compared  to  that  with  split  particles.  The 
effect  of  spatial  distributions  on  damage  evolution  is  studied  through  marked  correlation  functions 
M{r)  introduced  in  Pyrz  [22]  and  used  in  [27,  28].  Two  marks  associated  with  each  particle  are 
considered  for  their  relevance  to  damage  evolution.  They  are:  (a)  a  parameter  Rpg  which  is  defined 
as  the  ratio  of  the  maximum  principal  stress  to  the  critical  failure  stress  acr  for  an  undamaged 
particle,  and  as  the  ratio  of  the  current  overall  strain  to  the  strain  at  which  the  particle  had 
cracked,  for  a  cracked  particle  (^=  1).  Rps  signifies  the  propensity  to  advance  the  damage  state 
in  the  microstructure;  and  (b)  the  average  effective  plastic  strain  in  each  Voronoi  cell,  which 
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Second  order  intensity  function  K(r)  Cumulative  distribution  function  Cumulative  Distribution  Function 


(a) 


(e) 


(b) 


(f) 


Figure  3.4:  (a)  Cumulative  distribution  and  (b)  Probablity  density  functions  of  local  area  fractions; 
(c)  Cumulative  distribution  and  (d)  Probablity  density  functions  of  nearest  neighbor  distance;  (e) 
Second  order  intensity  and  (f)  Pair  distribution  functions  for  computer  simulated  microstructure. 
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Figure  3.5:  Damage  configurations  of  hard  core  microstructure  at  (a)  Ixx  =  0.8%,  (b)  —  1.4% 

and  (c)  Zxx  =  2.0%. Damage  configurations  of  triple  cluster  microstructure  at  (d)  Zxx  =  0.8%, 
(e)  Zxx  =  1.4%  and  (f)  Cxx  =  2.0%;  Corresponding  effective  plastic  strains(%)  for  triple  cluster 
microstucture  at  (g)  Zxx  =  1-4%  and  (h)  Zxx  =  2i^%,  for  particle  splitting  with  (Tct  =  300MPa. 


300.0 


Macroscopic  Tensile  Strain% 

Figure  3.6:  Macroscopic  stress-strain  response  for  computer  generated  microstructures  due  to  par¬ 
ticle  cracking  and  splitting. 

characterizes  evolving  matrix  failure  due  to  presence  of  damaged  particles. 

The  marked  correlation  function  M{r),  for  which  a  mathematical  formula  is  presented  in  [22], 
establishes  the  effect  of  microstructural  morphology  on  the  mark.  M{r)  for  the  two  patterns  and 
for  a  uniform  microstructure  (M(r)=l)  are  compared  in  figure  3.4.1.  For  the  hard  core  pattern  the 
functions  rapidly  decay  to  unity,  but  for  the  clustered  patterns  the  decay  is  considerably  slower. 
Higher  M(r)  values  for  clustered  patterns  at  short  range  sampling  distances  r  represent  larger 
influence  of  the  damage  marks.  The  high  value  of  M{r)  for  cp  at  short  sampling  distances  indicates 
severe  matrix  straining  near  damaged  particles.  The  slower  attenuation  of  M(r)  for  Rpg  at  short 
to  medium  range  indicates  that  particle  cracking  is  a  major  mode  of  damage  evolution.  This 
function  is  effective  in  understanding  the  sensitivity  of  damage  variables  to  local  perturbations  in 
the  morphology,  and  can  provide  a  criterion  for  determining  the  optimal  RME  size. 

3.4.2  Particle  splitting  simulation  with  actual  micrographs 

In  this  example,  VCFEM  analyses  is  conducted  with  micrographs  obtained  from  serial  sectioning 
of  reinforced  Al-Si-Mg  alloys  containing  w  10%  or  20%  by  volume  of  Si  particulates  (see  [29]). 
The  material  is  developed  at  ALCOA  Technical  Center  by  rapid  solidification  of  fine  powders  using 
a  gas  atomization  process  [22],  to  achieve  equiaxed  Si  particles.  The  powder  is  consolidated  by 
cold  isostatic  compaction,  canned  and  degassed  at  454°C,  and  finally  consolidated  to  full  density 
by  hot  isostatic  pressing.  Two  types  of  microstructure  are  considered,  viz.  (a)  a  naturally  aged 
20%  volume  fraction  composite  with  mean  Si  particle  size  of  4.4  /im,  and  (b)  a  naturally  aged  10% 
volume  fraction  composite  with  mean  Si  particle  size  of  3.9  /xm.  Serial-sectioning  of  the  specimens 
yield  a  series  of  2-D  sections  as  discussed  in  [29],  which  are  then  digitized. 
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Marked  Correlation  Function 


(a)  (b) 

Figure  3.7;  Marked  correlation  function  for  (a)  ratio  of  maximum  principal  stress  to  critical  fail¬ 
ure  stress  {o’er)  (b)  average  effective  plastic  strain  in  each  Voronoi  cells,  for  computer  simulated 
microstructure. 
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Characterization 

Equivalent  microstructures  that  closely  approximate  the  actual  2-D  morphology  of  micrographs  and 
yet  are  computationally  less  intensive  are  generated.  In  this  process  each  particle  of  arbitrary  shape 
is  replaced  by  an  equivalent  ellipse,  constructed  by  equating  the  0-th,  1-st  and  2-nd  moments  of  the 
actual  particles  with  those  of  the  equivalent  ellipses.  The  moments  of  actual  particle  are  computed 
as  the  sum  of  moments  of  each  pixel  contained  within  the  particle,  while  the  moments  of  the  ellipse 
are  represented  in  terms  of  standard  geometrical  parameters.  The  procedure  yields  (i)  the  centroid 
(^c5?/c)»  (ii)  lengths  (a,  6)  of  the  major  and  minor  axes,  and  (iii)  angular  orientation  9  of  the  major 
axis  of  the  equivalent  ellipse,  details  of  which  are  discussed  in  [21].  An  optical  micrograph  of  a 
section,  overlapped  simulated  and  exact  microstructures,  and  the  Voronoi  cell  mesh  obtained  by 
surface  to  surface  tessellation,  are  presented  in  figures  (3.4.2a,  b,  and  c)  respectively. 

Two  sections  of  each  volume  fraction  10%  and  20%  are  analyzed.  The  Vf  =  10%  sections  have 
77  and  89  Si  particles,  while  the  V/  =  20%  sections  contain  97  and  106  Si  particles.  A  majority  of 
computer  results  are  explained  with  respect  to  sections  of  the  20%Vf,  for  which  the  microstructural 
element  has  dimensions  of  205/i  x  180^,  and  the  two  dimensional  area  fractions  for  the  sections 
are  calculated  to  be  »  A/  =  18.6%  and  A/  =  18.9%  respectively  (note  that  the  3-D  V/  ~  20%). 
Particle  size  distribution  histograms  (not  shown)  show  considerable  scatter  within  each  section  and 
also  between  sections.  Cumulative  (F)  and  probability  density  distribution  (/)  functions  of  the 
local  area  fraction  (A)  and  near-neighbor  distance  (d)  are  plotted  in  figures  3.4.1(a,b,c  and  d).  A 
comparison  of  the  distribution  functions  in  figures  3.4.1  reveal  that  the  particle  distribution  is  more 
in  line  with  the  hard-core  pattern.  Similar  observations  are  also  made  in  the  plots  of  the  second 
order  intensity  function  F(r)  and  the  pair  distribution  function  in  figure  3.4.1,  where  the  absence 
of  local  peaks  in  g(r)  signals  negligible  clustering.  These  observation  is  consistent  with  the  material 
fabrication  process  innwhich  the  Si  particulates  are  randomly  precipitated  from  the  mixture. 

Damage  Simulation 

Material  properties  of  the  constituents  are:  For  the  Al-Si-Mg  matrix:  Young’s  modulus  E  —  69 
GPa,  Poisson’s  ratio  u  =  0.33,  and  the  post  yield  elastic-plastic  behavior  (non-linear  isotropic 
hardening)  is  obtained  from  data  on  T4-matrix  presented  in  figure  (8)  of  Kiser  et.al.  [24].  For 
the  Si  inclusions:  Young’s  modulus  E  =  161  GPa,  Poisson’s  ratio  i/  =  0.2  and  mode  I  critical 
stress  intensity  factor  for  pure  Si  is  found  to  be  Kjc  —  0.6  MPa-y/m.  In  the  fracture  mechanics 
approach  for  determining  the  size  dependent  critical  fracture  stress  <7^  in  equation  (3.9),  the  initial 
flaw  size  c  is  assumed  to  be  proportional  to  the  average  equivalent  particle  diameter  of  Davg-  The 
proportionality  constant  e  is  calibrated  by  analyses  of  auxiliary  RMEs  created  from  micrographs 
of  other  sections  of  the  specimen.  A  comparison  is  made  between  the  computer  simulations  and 
experimental  observations  with  micrographs  for  (a)  the  number  of  cracked  particles  and  (b)  overall 
stress-strain  behavior.  The  estimate  is  obtained  to  be  e  =  p  =  .125  or  12.5%,  and  therefore  the 
critical  stress  to  fracture  is  taken  to  be  <Tct  =  ^  ■  For  the  approach  with  WeibuU  distribution, 

the  two  micrographs  of  10%V/  are  used  for  calculating  parameters  Ocr  and  m,  since  they  exhibit 
the  onset  of  particle  cracking  at  =  3%.  The  cross-sectional  area  A  of  each  individual  particle 
is  calculated  from  the  data  on  ellipses  in  figure  3.4.2(d).  The  maximum  principal  stress  aj  at 
(xx  =  3%  is  obtained  from  V CFE  analysis  without  any  particle  damage.  The  probability  of  failure 
P/(A,<t)  in  equation  (3.10)  is  assumed  to  be  >  0.95.  The  WeibuU  parameters  are  evaluated  to  be 
m  =  2.37  and  Oct  =  2.12  GPa,  by  comparison  with  the  micrograph  observations  and  by  fitting  the 
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Figure  3.8:  (a)  Optical  micrograph  of  a  section  of  Al-Mg-Si  composite  (Si  Vf  =  20%)  (b)  simulated 
microstructure  superimposed  on  the  micrograph  (c)  Voronoi  cell  mesh  resulting  from  tessellation. 


data  in  a  least-square  sense.  The  probability  function  Pf  of  individual  particles  indicate  that  both 
particle  size  and  stress  levels  contribute  independently  towards  cracking.  A  qualitative  comparison 
between  the  two  approaches  for  evaluating  critical  stress,  is  made  by  VCFEM  analysis  with  the 
20% V/  section  of  97  particles  strained  to  €xx  =  6%.  Progression  of  damage  in  the  microstructure 
with  increasing  strain  is  shown  in  figure  (3.4.2a, b  and  c).  Particles  colored  black  are  cracked  while 
the  grey  particles  are  uncracked.  A  comparison  of  the  fraction  of  different  sized  particles  that 
are  cracked  by  each  approach,  is  made  with  actual  micrographic  observation  (figure  3.4.2a)  in 
histograms  (3.4.2d  and  e).  The  fractions  for  the  actual  micrograph  are  shown  in  grey  with  dashed 
outline.  The  histograms  indicate  that  a  while  both  criteria  are  good  for  large  particles,  the  WeibuU 
distribution  based  approach  provides  a  better  agreement  with  the  micrographs  at  the  smaller  size 
range.  Hence,  it  is  used  in  aU  subsequent  simulations. 

A  damage  parameter  is  defined  as  /)  =  %ea  'oj‘Ml  ^PartfdL ’  ^^ich  accounts  for  size,  is  plot¬ 
ted  in  figure  3.4.2a  for  each  of  the  two  sections  of  the  F/  =  10%  &  20%  composites  as  a  function 
of  increasing  strain.  Experimental  values  of  the  corresponding  damage  parameter  are  given  in 
Hunt[22],  where  the  arccis  are  calculated  by  sectioning  after  straining  to  a  certain  level.  Conse¬ 
quently  a  single  data  point  is  obtained  from  each  specimen.  Three  experimental  data  points  are 
plotted  for  the  Vp  =  20%  composite,  while  the  Vp  =  10%  has  a  single  data  point.  Generally 
a  good  agreement  is  noted  between  simulated  results  and  experimental  data  for  the  Vj  —  20% 
composite.  Stress-strain  response  of  nearly  undamaged  microstructures  are  shown  in  figure  3.4.2b 
by  loading  RMEs  in  compression  to  an  average  strain  of  c^x  —  6%.  While  the  elastic  response 
are  not  very  different,  the  20%^/  composite  has  higher  yield  stress  and  higher  flow  stress  than  the 
10%F/  composite.  The  tensile  response  with  particles  cracking  according  the  WeibuU  criterion  is 
shown  in  figure  3.4.2  c.  Results  of  VCFEM  analyses  for  each  volume  fraction  are  averaged  over  the 
two  sections  and  are  compared  with  experimental  results  of  Kiser  and  Zok  [24].  WhUe  VCFEM 
analysis  is  in  2-D  and  the  experimental  results  are  for  3-D,  the  comparisons  have  a  good  quaUtative 
agreement.  The  cross-over  point  at  which  the  20%Vf  composite  becomes  less  stronger  than  the 
10%^/  composite,  is  approximately  at  €xx  ^  1-2  —  1.8%,  and  compares  weU  with  the  experimental 
value  of  Cxx  ^  1.8%.  The  stress  capacity  is  in  general  higher  for  VCFEM  predictions,  which  may 
be  attributed  to  the  constrained  plastic  flow  arising  from  plane  strain  constraints.  AdditionaUy, 
the  present  simulation  does  not  aUow  matrix  softening  which  can  also  lower  the  load  capacity.  The 
stress-strain  behavior  of  a  uniform  (square  edge)  microstructure  with  a  single  circular  inclusion 
of  volume  fractions  10%  and  20%  are  also  plotted.  The  10%  uniform  composite  does  not  crack 
for  the  range  of  strains  considered  and  predicts  a  stiff  response.  Failure  of  the  single  particle  in 
the  20%Fy  composite  results  in  an  abrupt  drop  in  load  capacity  and  yields  unreasonable  predictions. 

Contour  plots  of  particle  failure  probability  and  effective  plastic  strains  for  the  20%Vf  composite 
are  iUustrated  in  figures  (3.4.2)  and  (3.4.2)  respectively.  Damaged  particles  are  in  white  with  cracks 
in  figure  (3.4.2),  and  the  contour  plots  are  for  undamaged  particles  indicating  the  likelihood  of 
damage.  An  interesting  observation  made  from  these  plots  is  that  some  large  particles  which  exhibit 
a  higher  tendency  to  crack  at  early  stages  of  loading,  may  remain  intact  throughout  due  to  failure  of 
neighboring  particles  and  load  redistribution.  This  phenomenon,  also  noticed  with  particle  clusters 
in  the  previous  example,  illustrates  the  influence  of  evolving  microstructural  morphology  on  the 
propagation  of  damage.  A  particle  crack  induces  large  plastic  flow  in  the  neighboring  matrix  which 
causes  the  stress  to  rise  in  particles  in  this  region  and  eventually  initiate  a  crack.  The  plastic  strain 
distribution  in  figure  (3.4.2  a)  shows  localized  bands  of  severe  deformations  emanating  from  crack 
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Figure  3.9:  Simulated  configurations  of  evolving  damage  at  (a)  e^x  =  3.0%,  (b)  e^x  =  4.2%  and 
(c)  ixx  =  6.0%  with  a  WeibuU  distribution  based  damage  criterion;  (d,e)  Histograms  of  number 
of  damaged  particles  at  €yy  =  6%,  by  WeibuU  distribution  and  fracture  mechanics  based  damage 
criteria  respectively. 


24 


(a) 


1 

2 


400.0 


300.0 


0^ 

g 

I 

C/5 

O 


2  200.0 


H 

o 

o 


100.0 


0.0 


Matrix 

10%  Composite  (Kiser  et.al.  [26]) 
20%  Composite  (Kiser  et.al.  [26]) 
10%  Composite  (VCFEM  averaged) 
20%  Composite  (VCFEM  averaged) 
10%  Uniform  Composite  (VCFEM) 
20%  Uniform  Composite  (VCFEM) 


0.0  1.0  2.0  3.0  4.0  5.0 

Macroscopic  Tensile  Strain% 


6.0 


7.0 


(b)  (c) 

Figure  3.10:  (a)  Area  fraction  of  cracked  particles  as  a  function  of  the  macroscopic  strain.  Macro¬ 
scopic  stress-strain  response  for  the  Al-Si-Mg  composite  microstructures  at  two  volume  fractions 
(b)  compressive  response  (c)  tensile  response. 


25 


Max. 


qn  V 

cHE^  On 

g^»®» 

(p  %  ^ 

rnN  m  (po 


# 

#  • 


•Cb 

^0) 


■Q) 


0^ 


#  f%' 


(D 


(j^  _  % 


QX  Ik 

^(6  ®(P.  V. 


m 


®®rc^  -cd'^^'  »- 

cf  2p  ^ 

»®  ^  %  .  %  Q3»> 

*  Q)®..*Cb®®.  ' 

..  CP  •c^,»_* 


(a)  (b) 

Figure  3.11:  Contour  plots  of  particle  fracture  probability  by  the  Weibull  damage  criterion  im¬ 
mediately  before  (a)  Ixx  =  4.2%  and  (b)  I^x  =  6.0%  for  the  Al-Si-Mg  composite  microstructure 
{Vf  =  20%).  (Damaged  particles  are  in  white  with  cracks). 


tips  and  propagating  to  neighboring  particles  with  cracks.  The  remainder  of  the  matrix  undergoes 
relatively  smaller  deformations.  Marked  correlation  functions  M{r)  are  plotted  in  figure  (3.4.2)  as 
functions  of  distance  and  particle  shapes,  for  the  two  sections  of  20%V/  composite.  Two  marks 
viz.  (a)  particle  fracture  probability  Pf  and  (b)  effective  plastic  strain  ^  are  selected  as  qualitative 
indicators  of  microstructural  damage.  Uniform  M{r)  plots  of  unit  value  occur  for  regular  patterns 
and  correspond  to  identical  marks.  Figure  (3.4.2a  and  b)  of  M(r)  with  respect  to  the  sampling 
distance  r  shows  that  the  functions  for  Pf  quickly  stabilize  near  the  unit  value,  while  the  decay 
is  slower  for  €?  with  a  few  abrupt  peaks.  The  lack  of  strong  clustering  in  these  patterns  leads  to 
smaller  influence  regions  for  these  microstructures.  Marked  correlation  functions  with  respect  to 
the  relative  difference  in  form  factors  Ff  are  plotted  in  figure  (3.4.2c  and  d).  The  form  factor  is  an 
indicator  of  the  deviation  in  shape  from  a  perfect  circle  (F/=l),  and  for  an  elliptical  shape,  may 
be  expressed  as  (see  [28]): 


Ff  =  — — 

perimeter"^ 


R  =  y/^  ;  perimeter  «  7r[1.5(a  +b)-  V^]  (3.11) 


where  a  and  b  are  the  major  and  minor  axes.  The  maximum  observed  form  factor  Ff  is  computed 
to  be  0.97  while  the  minimum  is  0.48.  The  figures  (3.4.2c  and  d)  depict  increasing  correlation 
functions  especially  for  This  infers  that  shape  has  significant  influence  on  the  damage  in  these 
microstructures. 
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Figure  3.12:  EflFective  plastic  strain  (%)  contours  at  (a)  Ixx  =  4.2%  and  (b)  =  6.0%  for  the 

Al-Si-Mg  composite  microstructure  (V/  =  20%). 

3.5  Conclusions 

This  work  is  devoted  to  stress  analysis  of  non-uniform,  ductile  matrix  composite  microstructures 
with  particle  cracking,  by  the  Voronoi  CeU  finite  element  method  (VCFEM).  The  computational 
model  assumes  complete  cracking  at  the  onset  of  damage,  and  differentiates  between  the  behavior 
of  fuUy  cracked  particles  and  split  particles.  The  uniqueness  of  this  model  lies  in  its  ability  to 
model,  continuously  changing  element  topology  due  to  progressive  material  failure,  with  no  user 
interference.  Validation  of  the  computation  model  for  damage  is  done  through  various  studies, 
including  comparison  with  other  numerical  studies  in  the  literature  that  use  conventional  finite 
element  codes.  These  studies  have  predominantly  analyzed  simple  uniform  distributions  with  pre¬ 
existing  cracks.  Good  agreement  is  obtained  in  these  comparison  studies,  both  from  a  macroscopic 
and  microscopic  point  of  view. 

A  major  advantage  of  VCFEM  is  that  it  can  be  used  for  analyzing  damage  in  nonuniform  real 
micrographs  without  making  major  morphological  simplifications.  A  set  of  computer  generated 
hard-core  and  clustered  microstructures  are  simulated  to  understand  the  effect  of  spatial  distribu¬ 
tion  on  damage  evolution.  Damage  initiates  within  each  cluster,  propagates  within  the  cluster  and 
finally  links  up  with  damage  in  the  neighboring  clusters.  Regions  of  severe  plastic  flow  exist  in 
the  matrix  ahead  of  split  particles,  indicating  possible  sites  of  matrix  failure.  Particle  splitting  is 
found  to  yield  much  softer  overall  response  than  particle  cracking.  In  a  concluding  example,  the 
VCFEM  model  is  directly  constructed  from  a  digitized  optical  micrograph  of  an  Al-Si-Mg  compos- 
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Figure  3.13:  Marked  correlation  functions  of  (a)  failure  probability  and  (b)  average  effective  plastic 
strain  in  each  Voronoi  ceU  as  a  function  of  radial  distance,  and  marked  correlation  functions  of  (c) 
failure  probability  and  (d)  average  effective  plastic  strain  in  each  Voronoi  ceU  as  a  function  of  form 
factor;  for  the  Al-Si-Mg  composite  microstructure  (V/  =  20%). 
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ite  system.  VCFEM  results  are  compared  with  results  of  experiments  with  this  material  for  two 
volume  fractions.  Good  concurrence  is  obtained  in  the  overall  stress-strain  behavior,  as  well  as  in 
the  number  of  damaged  particles  in  the  microstructure.  Larger  particles  tend  to  fail  at  lower  load 
levels  and  therefore  serve  as  sites  of  damage  initiation.  Nevertheless,  smaller  particles  may  also 
damage  with  loading,  due  to  stereological  factors  like  proximity  with  other  particles  and  relative 
shape.  It  is  noted  that  damage  evolution  in  real  microstructures  is  a  gradual  process,  and  takes 
place  by  progressive  particle  failure.  Material  behavior  is  therefore  misrepresented  with  the  single 
cell  models  where  particle  cracking  result  in  abrupt  changes  in  response.  Statistical  descriptors  are 
used  as  indicators  of  morphological  influence  on  the  damage  state.  The  efficiency  of  the  VCFEM 
codes  is  noteworthy.  A  comparison  with  conventional  FEM  packages  for  undamaged-nonuniform 
and  damaged-uniform  materials  show  a  ~  30-50  times  reduction  in  computing  time.  Currently, 
matrix  cracking  phenomenon  is  being  incorporated  in  VCFEM  and  this  will  be  reported  in  future. 
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Chapter  4 


Multi-level  Computational  Model 
For  Multiscale  Damage  Analysis  In 
Composite  And  Porous  Materials 

4.1  Introduction 

Heterogeneous  structures  with  second  phase  inclusions  or  voids  in  the  microstructure  are  conven¬ 
tionally  analyzed  with  macroscopic  properties  obtained  from  homogenization  of  response  at  smaller 
(meso-,  micro-)  length  scales.  The  mathematical  homogenization  theory,  which  uses  asymptotic 
expansions  of  displacement,  strain  and  stress  fields  about  macroscopic  values,  has  been  used  as  a 
tool  for  analyzing  multiple  scale  responses  in  [Ij  2,  3,  4].  The  method  is  based  on  eissumptions  of 
spatial  periodicity  of  microscopic  representative  volume  elements  (RVE)  and  local  uniformity  of 
macroscopic  fields  within  each  RVE.  It  decomposes  the  multiscale  boundary  value  problem  into  a 
decoupled  set  of  micro-scale  RVE  problem  and  a  macro-scale  problem.  Concurrent  finite  element 
analyses  are  executed  at  the  each  scale  for  information  transfer  between  the  scales.  Multiple  scale 
analysis  of  linear  elastic  reinforced  composites  by  this  method  have  been  conducted  by  Fish  et. 
al.  [5,  6],  Kikuchi  et.  al.  [7,  8].  For  nonlinear  materials,  the  homogenization  methods  have  been 
extended  by  Suquet  [9],  Fish  et.al.  [10],  Guedes  [11]  and  Cheng  [12].  The  method  has  also  been 
implemented  to  simulate  damage  by  fiber-matrix  debonding  in  linear  elastostatics  [13]  and  fiber 
rupture  using  a  phenomenological  damage  model  [14]. 

Despite  its  advantages,  asymptotic  homogenization  has  suffered  shortcomings  arising  from  effi¬ 
ciency  and  accuracy  considerations.  Enormous  computational  efforts  can  result  with  this  method 
due  to  the  fact  that  at  each  integration  point  in  the  macroscopic  model,  boundary  value  problems 
of  the  microstructural  RVE  should  be  solved  twice.  To  economize  computations,  many  studies  have 
assumed  simple  unit  cells  models  of  the  microscopic  RVE.  Such  idealizations  may  however  be  unre¬ 
alistic  for  deformation  and  failure  analysis  of  many  materials.  The  homogenization  method  has  an¬ 
other  major  limitation  stemming  from  its  basic  assumptions,  viz.  (a)  uniformity  of  the  macroscopic 
fields  within  each  RVE  and  (b)  spatial  periodicity  of  the  RVE.  The  uniformity  assumption  is  not 
appropriate  in  critical  regions  of  high  gradients,  where  the  macroscopic  fields  can  vary  considerably. 
Free  edges,  interfaces,  macrocracks,  neighborhood  of  material  discontinuities  and  most  importantly 
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in  the  regions  of  evolving  microscopic  damage  and  instability  are  potential  sites  of  nonuniformity. 
Furthermore,  statistical  periodicity  implies  that  the  RVE  may  be  repeated  to  represent  the  entire 
neighborhood  of  a  macroscopic  point.  For  non-uniform  microstructures  sufficiently  large  portions 
should  be  considered  as  RVE  for  homogenization  analysis.  Unit  cell  models  are  severely  limited  in 
this  respect.  Even  higher  order  theories  of  homogenization  may  be  computationally  un viable.  A 
few  effective  global  local  techniques  based  on  hierarchical  decomposition  and  superposition  of  field 
variables  have  been  proposed  by  Belytschko  [15],  Reddy  [16]  and  Hughes  [17].  Pagano  and  Rybicki 
[19]  had  discussed  the  breakdown  of  effective  modulus  theory  for  composite  laminates  with  free 
edges  and  the  need  for  global-local  techniques.  Fish  et.al.  [6,  18]  have  used  global-local  techniques 
with  multigrid  methods  to  extend  the  multiple  scale  modeling  to  non-periodic  materials.  Zohdi 
and  Oden  [20,  21]  have  developed  a  homogenized  Dirichlet  projection  method  (HDPM)  which  re¬ 
solves  the  microstructural  effects  at  different  scales  on  the  macroscopic  response  of  heterogeneous 
structures. 

To  improve  the  accuracy  of  multiple  scale  analyses  involving  microstructural  damage  of  elastic- 
plastic  composite  and  porous  structures  without  computational  efficiency,  an  adaptive  multi-level 
method  is  proposed  in  this  work.  It  uses  computational  hierarchy  in  the  different  levels  to  con¬ 
currently  predict  evolution  of  variables  at  the  structural  and  microstructural  scales,  as  weU  as  to 
track  the  incidence  and  propagation  of  damage.  Analysis  of  microstructural  response  with  arbi¬ 
trary  distributions,  shapes  and  sizes  of  heterogeneities  is  conveniently  done  by  the  Voronoi  Cell 
finite  element  model  (VCFEM)  [22,  23,  24,  25].  The  VCFE  model  naturally  evolves  by  tessellation 
of  the  microstructure,  to  generate  a  morphology  based  network  of  multi-sided  Voronoi  polygons. 
Each  Voronoi  ceU  with  the  embedded  heterogeneity  is  treated  as  a  FEM  element  in  this  model. 
Incorporation  of  micromechanics  based  assumptions  into  a  hybrid  finite  element  formulation  im¬ 
parts  a  high  level  of  computational  efficiency  with  sufficient  accuracy  and  resolution  in  this  method. 
Furthermore,  pre-processing  efforts  in  generating  microstructural  models  are  drastically  reduced. 
Various  thermo-elastic  and  elastic-plastic  problems  of  composite  and  porous  materials  have  been 
successfully  analyzed  by  VCFEM  [22,  24,  25].  Progressive  damage  by  particle  cracking  has  been 
done  in  [23]  and  VCFEM  has  been  used  to  relate  image  analysis  and  quantitative  characterization 
with  microstructural  response  in  [24].  For  periodic  representative  volume  elements  (RVE)  of  elastic 
and  elastic-plastic  materials,  the  microstructural  VCFEM  has  been  coupled  with  structural  analysis 
codes  by  using  asymptotic  homogenization  in  [26,  27].  When  the  method  fails  due  to  questionable 
eissumptions  of  macroscopic  uniformity  and  statistical  periodicity,  a  combination  of  homogeniza¬ 
tion  and  global-local  methods  is  a  necessity.  This  implementation  is  a  nontrivial  undertaking  due 
to  lack  of  apriori  knowledge  of  regions  requiring  differential  resolution.  The  adaptive  multi-level 
methodology  developed  in  this  work  adresses  this  challenge.  The  model  differentiates  between 
non-critical  and  critical  regions  for  which  the  model  ranges  from  macroscopic  computations  using 
continuum  constitutive  relations  to  zooming  in  at  ‘hotspots’  for  pure  microscopic  simulations.  The 
adaptive  process  significantly  increases  the  efficiency  while  retaining  appropriate  level  of  accuracy 
in  the  hierarchy.  The  chapter  begins  with  a  discussion  on  two  scale  computations.  It  introduces 
three  levels  of  the  computational  domain  which  include  macro,  macro-micro  and  microscopic  anal¬ 
ysis.  Numerical  examples  are  conducted  for  heterogeneous  materials  with  a  variety  of  microscopic 
architectures  to  support  its  development. 
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Figure  4.1:  A  heterogeneous  structure  showing  various  scales:  (a)  The  structure  at  the  macroscopic 
scale  of  applied  loads  (b)  A  representative  volume  element  (RVE)  at  the  microscopic  scale  with 
the  VCFE  model  and  (c)  A  Voronoi  cell  element  at  the  scale  of  a  single  heterogeneity  or  basic 
structural  element. 
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4.2  Two  Way  Coupling  for  Multiple  Scale  Modeling 

Multiple  scale  modeling  of  heterogeneous  materials  is  necessary  to  concurrently  account  for  spatial 
variability  at  the  macro-  and  micro-scales.  An  effective  model  of  this  class  requires  two-way  coupling 
for  efficient  computing,  as  well  as  for  accurate  representation  of  the  necessary  variables  at  different 
scales.  The  first  is  a  ‘bottom  up’  coupling  for  determination  of  equivalent  homogeneous  behavior  at  a 
macroscopic  point  x,  as  a  function  of  the  microstructural  geometry  and  behavior  of  the  constitutive 
phases,  but  independent  of  applied  loads  to  the  structure.  In  the  homogenization  procedure,  an 
isolated  representative  volume  element  (RVE)  Y(x)  C  3?^  is  identified  at  microstructural  scale  of 
heterogeneities  (figure  4.1b).  The  scale  y  of  the  RVE  domain  Y(x)  may  be  large  with  respect  to 
the  characteristic  length  1  of  microscopic  heterogeneities,  but  is  significantly  small  compared  to  the 
the  macroscopic  length  scale  L  of  the  structure  and  applied  loads.  Homogenized  variables  at  the 
macroscopic  scale  are  obtained  by  volume  averaging  of  variables  in  the  RVE,  following  the  definition 

<  f  >Y=  ^  J^f{y)dy  (4.1) 

The  condition  for  macroscopic  homogeneity,  according  to  the  Hill-Mandel  hypothesis  [29],  assumes 
equivalence  of  strain  energy  for  the  actual  and  equivalent  homogenized  media.  Thus  for  a  statically 
admissible  stress  field  <r(y)  and  kinematically  admissible  strain  field  e(y), 

<  <r  :  €  >Y=<  o- >y:<  e>Y  Vy  €  Y  (4.2) 

The  microscopic  stress  <r(y)  and  strain  e(y)  fields  satisfying  the  homogeneity  condition  (4.2)  may 
be  obtained  by  solving  boundary  value  problems  for  the  RVE  Y  with  prescribed  homogeneous 
stress/strain  or  periodicity  boundary  conditions,  stated  as: 

T*"  =<  or  >y  •n(y)  =  o-  •  n(y)  on  dY  :  Uniform  Traction  (a) 

u*’  =<  6  >Y  -y  on  dY  :  Uniform  Strain  (b) 

u**  =<  €  >Y  -y  +  u  =<  6  >Y  -y  +  u  +  kY  on  dY  :  Y-Periodicity  (c)  (4-3) 

where  k  is  an  integer  and  Y  is  the  basic  period  of  the  Y-periodic  displacement  functions.  The 
macroscopic  constitutive  equations  are  obtained  by  solving  a  boundary  value  problem  of  the  RVE 
Y  with  one  of  the  three  sets  of  boundary  conditions  in  equation  (4.3),  followed  by  the  averaging 
process  in  equation  (4.1).  For  linear  elastic  constituent  phases  in  Y,  the  relation  between  the  strain 
energy  functions  has  been  established  in  Suquet  [30]  as: 

<  e  >:  E*  :<  e  >  <  <  e  >:  eJ,,  :<  e  >  <  ,<  e  >:  E*,  :<  e  >  VE  C  (4.4) 

where  E^,.,  Ep^,.,  Ej^j^  ^^e  respectively  the  homogenized  stiffness  tensors  evaluated  with  uniform  trac¬ 
tion,  periodicity  and  uniform  strain  boundary  conditions,  and  <  e  >  is  the  macroscopic  (applied  or 
averaged)  strain  field.  The  discrepancy  with  the  kinematic  and  kinetic  boundary  conditions  reduce 
with  increasing  size  of  Y.  It  is  generally  concluded  [8],  that  for  the  same  RVE  size,  the  periodicity 
boundary  conditions  are  expected  to  yield  more  accurate  statistically  homogenized  constitutive 
parameters  and  macroscopic  properties. 
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The  other  coupling  is  the  ‘top  down’  where  the  evolution  of  variables  are  evaluated  in  the 
microstructure  from  known  macroscopic  variables,  by  a  process  termed  as  localization.  In  those 
regions,  where  the  microstructure  may  admit  a  RVE  Y,  the  microscopic  variables  can  be  evaluated 
by  solving  a  boundary  value  problem  with  imposed  macroscopic  strains  and  the  local  periodicity 
condition  in  equation  4.3c.  In  other  regions,  where  assumptions  of  local  periodicity  of  the  RVE 
may  be  unrealistic,  the  localization  process  will  entail  direct  interfacing  of  the  microstructural  and 
macroscopic  regions. 

4.3  A  Multi  Level  Model  for  Coupling  Different  Scales 

In  the  spirit  of  true  two  way  coupling  of  multiple  scale  problems,  the  computational  domain  is 
adaptively  decomposed  into  three  levels  of  hierarchy  based  on  requirements  of  resolution.  Such 
hierarchy  is  intended  to  increase  computational  efficiency  as  weU  as  accuracy  in  concurrent  predic¬ 
tion  of  variables  at  the  continuum  and  microstructural  scales.  As  proposed  in  [31],  the  model  uses 
homogenization  of  microstructural  RVE  solutions  to  evaluate  homogenized  properties  and  cascades 
down  the  scales  at  hotspots  of  evolving  damage.  The  three  levels  of  hierarchy  with  requirements  of 
increasing  resolution  (figure  4.2)  are  as  follows. 

•i.  Computational  Subdomain  Level-0; 

These  correspond  to  non-critical  macroscopic  regions  in  figure  4.2a,  where  deformation  variables 
are  relatively  uniform  and  periodicity  conditions  may  be  assumed  for  the  underlying  material  RVE. 
Scale  effects  are  negligible  in  this  region  and  local  constitutive  relations  may  be  derived  from  pos¬ 
tulates  of  the  RVE  approaching  zero  volume.  Continuum  level  anisotropic  plasticity  constitutive 
relations,  that  are  consistent  with  the  actual  microstructural  constitution,  are  developed  for  macro¬ 
scopic  modeling  of  these  regions. 

The  /eue/-0  macroscopic  simulations  are  accompanied  by  element  refinement  or  h-adaptation  for 
two  reasons.  The  first  is  to  identify  and  reduce  a  chosen  ‘error  measure’  in  the  macroscopic  com¬ 
putational  model.  A  second  attribute  is  that  it  enables  the  computational  model  to  ‘zoom  in’  on 
regions  of  evolving  nonuniformity  due  to  microscopic  non-homogeneity.  This  reduces  the  disparity 
in  size  of  the  macro-  and  micro-  scale  elements  by  successive  refinement  of  macroscopic  elements 
in  the  critical  regions  as  shown  in  figure  4.2a. 

•ii*  Computational  Subdomain  Level-1;  These  are  regions  that  face  imminent  microscopic 
non- homogeneity  and  resulting  macroscopic  nonuniformities  (figure  4.2a, b).  Though  the  computa¬ 
tions  are  still  macroscopic,  concurrent  monitoring  of  the  development  of  damage  and  instabilities  in 
the  RVE  is  possible  in  this  level.  For  concurrent  macro-  and  micro-  scale  analyses  the  asymptotic 
homogenization  methods,  which  is  bctsed  on  the  existence  of  an  RVE,  is  used.  Macroscopic  element 
refinement  by  h-adaptation  continues  for  this  level. 

•“i*  Computational  Subdomain  Level-2:  These  critical  regions  materialize  with  the  evolution 
of  microstructural  damage  in  the  form  of  evolved  microcracks  or  instabilities  (figure  4.2b),  leading  to 
high  macroscopic  field  gradients.  The  assumptions  of  macroscopic  uniformity  and  local  periodicity 
are  unrealistic.  To  realize  scale  effects,  it  is  required  that  the  level-1  macro-micro  computational 
model  switch  to  a  completely  microscopic  model,  encompassing  large  portions  of  the  microstructure. 
A  detailed  flow  chart  of  the  adaptive  hierarchical  process  is  depicted  in  figure  4.3. 
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Figure  4.2:  A  hierarchical  multi-level  computational  domain;  Level-0  for  macroscale  continuum 
modeling  (b)  Level- 1  for  coupled  macro-microscopic  (RVE)  modeling  with  asymptotic  homogeniza¬ 
tion  and  (c)  Level-2  for  pure  microscopic  modeling 
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The  substructured  computational  domain  is  delineated  as  an  elastic-plastic  body  of  material 
domain  ^mat  that  consists  of  regions  for  which  the  RVE  is  repeated  periodically,  and 

also  of  regions  ^MAT[np)  where  the  periodicity  assumptions  do  not  hold,  i.e. 

^MATi-X.,  u)  =  9.mAT(p){x,  u)  U  ^MAT(np){x,  u)  (4.5) 

The  macroscopic  regions  of  periodicity  are  further  constituted  of  repeating  a  large  number  of  RVE’s 

^MAT(p)  =  uS^Yf  where  vf ¥[  (4.6) 

Here  N{p)  corresponds  to  the  number  of  different  RVE’s  Y|.  in  the  periodic  regions,  and  for  aH 
practical  purposes,  oo  corresponds  to  a  sufficiently  large  number.  The  non-periodic  region  ^MAT(np) 
is  defined  as  the  set  of  aU  microstructural  regions  for  which  the  N{np)  RVE’s  are  not  repeated,  i.e. 

^MAT(np)  =  ufif  ^Yf  ;  Yfc  n  Y,  =  0  V  A:  ^  /  (4.7) 

The  level-0  and  level-1  of  the  computational  domain  ft  correspond  to  the  periodic  regions,  while 
the  level-2  belongs  to  the  non-periodic  regions  as 

QlO  U  fill  C  ^MAT{p)  ;  ^12  c  ^MAT{np)  (4.8) 

4.4  Homogenization  with  Voronoi  Cell  FEM 

4.4.1  Asymptotic  homogenization  with  microstructural  periodicity 

Consider  a  heterogeneous  structure  occupying  a  region  ^structure  (figure  4.1a),  for  which  the  hetero¬ 
geneous  microstructure  constitutes  of  spatially  repeated  RVE’s  Y(x)  about  a  macroscopic  point  x 
as  shown  in  figure  4.1b.  The  RVE  is  discretized  into  a  mesh  of  Voronoi  cells,  which  naturally  evolve 
from  the  microstructure  by  Dirichlet  tessellation.  In  the  Voronoi  cell  FEM  (VCFEM),  each  Voronoi 
cell  represents  a  basic  structural  element  (BSE)  which  is  the  neighborhood  of  a  heterogeneity  in 
the  microstructure.  The  dimensions  of  the  Y(x)  are  typically  very  small  in  comparison  with  the 
structural  dimensions  X,  i.e.  j  is  a  very  small  positive  number  e.  Due  to  variation  of  evolutionary 
variables  in  a  small  neighborhood  e  of  the  macroscopic  point  x,  all  variables  are  assumed  to  exhibit 
dependence  on  both  length  scales  i.e.  =  $(x,  -),  where  y  =  7.  The  superscript  e  denotes  associ¬ 
ation  of  the  function  with  the  two  length  scales  and  hence  corresponds  to  a  connected  structural 
and  microstructural  domain.  The  assumption  of  periodic  repetition  of  the  microstructure  about  x 
makes  the  dependence  of  the  function  on  y(=  7),  Y-periodic  [7,  4,  8,  14].  For  small  deformation 
elasto-plasticity,  the  rate  forms  of  the  equilibrium,  kinematic  and  constitutive  relations  are  given 
as: 

Equilibrium  :  =  — /j,  Kinematics  : 

z 

Constitutive  :  in  fi'  (4.9) 

where  &fj(x,y),  ifj(x,y)  and  «^(x,y)  are  Y-periodic  rates  of  stress,  strain  and  displacement  fields 
respectively.  Furthermore  the  boundary  tractions  and  displacements  are  assumed  to  satisfy  the 
following  prescribed  conditions. 

a-jUj  =  it  on  Ft  ,  w;  =  it;  on  Tt^  (4.10) 


(duk  ,  dul\ 
W  9x1 ) 
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where  n  is  the  unit  normal  to  the  boundary.  The  Y-periodic  displacement  rate  or  increment  field 
is  approximated  by  an  asymptotic  expansion  about  x  with  respect  to  the  parameter  e: 


=  «i(x,y)  +  ewi(x,y)  +  e'^uU^^y)  +  ' 


y  = 


(4.11) 


Noting  that  the  spatial  x'  derivative  of  any  function  depends  on  the  two  length  scales  and  is  given 


as: 


dx^-  V  ^  ’  €  )  dxi  €dyi 


the  stress  rate  tensor  can  be  expressed  as 


^  ^0  I  I  I  ^2  ^3 


(4.12) 


(4.13) 


where 


^0_pe  ^2  _  pe  (dul.dul\  .  . 

From  equations  4.9  and  4.13,  and  using  the  periodicity  condition  on  the  RVE  boundary  /  dTy  = 
0,  it  can  be  proved  (see  [26,  27])  that 


.5=0,  =  g  =  +  =  0 


(4.15) 


Furthermore,  by  neglecting  the  terms  associated  with  e  or  higher  in  equation  (4.13),  the  constitutive 
relation  in  the  Y-domain  is  expressed  as 


•£  -1  me  •£  r.£  \ 

^ij  —  ^ij  —  ^ijkl^kl  —  ^ijkli  +  Qy^  ) 


(4.16) 


Here  eh  is  the  microstructural  strain  rate  tensor,  for  which  is  an  averaged  macroscopic  part 

du^  .  . 

^  is  denoted  as  a  fluctuating  strain  rate  tensor  [9].  Due  to  linearity  of  the  rate  problem, 
uj  and  the  microscopic  equilibrium  condition  can  be  expressed  as 


=  K'(y^  _  0 


(4.17) 


In  equation  (4.17),  afj  is  a  Y-antiperiodic  function  and  is  a  Y-periodic  function  represent¬ 
ing  characteristic  modes  of  the  deformation  rate  in  the  RVE.  Substituting  equation  (4.17)  in  the 
constitutive  relations  of  equation  (4.9)  yields  the  microscopic  constitutive  relations  as: 


(4.18) 


where  Sij  is  Kronecker  delta.  The  mean  of  equation  (4.18)  yields  the  homogenized  elastic-plastic 
tangent  modulus  for  use  in  the  macroscopic  analysis,  in  the  form 

®{$r  =<  »i!  >y=  I  »fldY  =  lY]  X  ('‘•19) 
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(4.20) 


The  macroscopic  stress  and  strain  relation  can  thus  be  stated  as 

s«(x)  =<  >y=  Efj„juix) 

where  the  homogenized  variables  are  S{x)  =<  <r^(x,y)  >y  and  e(x)  =<  €^(x,y)  >y.  The 
incremental  small  deformation  analysis  for  elastic-plastic  materials  can  be  conducted  with  the 
homogenized  modulus  at  the  macroscopic  level  and  by  using  the  Voronoi  cell  finite  element  model 
(VCFEM)  for  solving  the  microscopic  problem. 


4.4.2  The  Voronoi  cell  FEM  for  microstructural  analysis 

The  Voronoi  ceU  finite  element  model  has  been  successfully  developed  for  composite  and  porous 
materials  in  [22,  23,  25].  Arbitrary  dispersion  patterns,  shapes  and  sizes  of  heterogeneities  are 
readily  modeled  by  VCFEM.  The  computational  model  naturally  evolves  by  Dirichlet  tessellation 
of  the  microstructure  as  shown  in  figure  4.1b.  Each  Voronoi  ceU  with  the  embedded  inclusion  or 
void  is  treated  as  an  element  in  this  formulation.  Preprocessing  efforts  are  drasticaUy  reduced,  as  a 
consequence.  In  [23,  28],  the  VCFEM  formulation  has  been  extended  to  include  damage  evolution 
in  the  form  of  paxticle  cracking,  where  the  crack  is  realized  as  an  eUiptical  void.  Each  Voronoi  ceU 
element  is  amenable  to  change  in  topology  from  two  constituent  phases  (matrix  and  inclusions)  in 
undamaged  ceUs,  to  three  phases  (matrix,  inclusion  and  crack)  in  damaged  ceUs.  Complete  particle 
cracking  or  spUtting  is  assumed  to  occur  at  the  very  onset  of  damage. 


The  Voronoi  ceU  finite  element  formrdation  constructs  a  hybrid  element  by  combining  the  aspects 
of  finite  element  methods  with  important  micromechanics  considerations.  Use  of  a  hybrid  stress 
based  formulation  results  in  a  high  level  of  accuracy  with  a  significantly  reduced  degree  of  freedom, 
compared  to  displacement  based  FEM  models.  Consider  a  typical  representative  volume  element  Y 
consisting  of  N  undamaged  and/or  damaged  particles,  that  are  contained  in  each  of  the  N  Voronoi 
ceU  elements,  as  shown  in  figure  4.1(b).  The  assumed  stress  hybrid  formulation  in  the  Voronoi 
ceU  finite  element  method  (VCFEM)  requires  independent  assumptions  of  an  equilibriated  stress 
field  (<r)  in  the  interior  of  each  element  Yg,  and  compatible  displacement  fields  u  on  the  element 
boundary  dYg,  u'  on  the  matrix-inclusion  interface  dYg  and  u"  on  the  crack  boundary  dYgr-  In  an 
incremental  formulation  for  elasto-plasticity,  the  incremental  variational  formulation  introduces  an 
element  energy  functional, 


Ilf(A<r,  Au)  =  —  /  AB(<r,Aa-)  dY  —  f  e  :  A<r  dY  -f- 

JYg  JYg 

I  {or  +  Act)  •  n®  •  (u  -h  Au)  ddY  (Inter-element  Traction  Reciprocity) 

JdYe 

~  (*  +  At)  •  (u  -f  Au)  dr  (Boundary  Traction) 

—  I  (o-”*  -f  Ao-^  -  <7®  —  Ao-®)  •  n®  •  (u'  -f  Au')  dY  (Matrix-Inclusion  Interface  Traction) 

—  I  (o’*^  +  A<t®)  •  n®'"  •  (u"  -f  Au")  dY  (Crack  Boundary  Traction)  (4-21) 

JdYcr 


where  AB  is  the  increment  of  complimentary  energy  density.  Variables  (<r,u)  correspond  to  values 
at  the  beginning  of  an  increment,  while  variables  (A<t,  Au)  are  the  corresponding  increments  in 
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a  load  increment  or  step.  Outward  normals  on  dYc  and  dYcr  are  denoted  by  n®,  and 
respectively.  Superscripts  m,  c  and  cr  are  associated  with  the  matrix,  inclusion  and  crack  phases 
respectively  in  each  Voronoi  cell  element.  The  total  energy  for  the  entire  RVE  of  N  Voronoi  cells 
is  obtained  as  11^  =  52^1  n^.  Setting  the  first  variation  of  11^  in  equation  (4.21)  with  respect  to 
stress  increments  Air  to  zero  yields  the  element  compatibility  as  the  Euler  equation,  while  setting 
the  first  variations  of  11^  with  respect  to  the  independent  boundary  displacements  Au,  Au'.  and 
Au"  to  zero,  yield  the  inter-element  boundary  traction  reciprocity,  matrix-inclusion  interface  trac¬ 
tion  reciprocity  and  zero  traction  crack  boundary  condition  respectively.  Independent  assumptions 
on  stress  increments  Acr  are  made  in  the  matrix  and  inclusion  phases  in  each  element,  thus  allowing 
stress  discontinuities  across  the  interface.  In  this  process  special  forms  of  the  Airy’s  stress  func¬ 
tion  ^(x,y)  to  enhance  computational  efficiency,  has  been  developed  in  [22,  23]  for  equilibriated 
stress  fields.  The  functions  facilitate  stress  concentration  near  the  interface  and  crack  boundary, 
accounting  for  the  shape  of  the  inclusion  and  crack  and  also  help  satisfy  traction  reciprocity  at  the 
interfaces  dY^  and  dYcr.  Furthermore,  they  decay  at  large  distances  from  the  interfaces.  Compati¬ 
ble  displacement  increments  are  generated  on  each  of  the  boundaries/interfaces  dYg,  dYc  and  dYcr 
by  interpolation  in  terms  of  generalized  nodal  displacements  as  , 

{Au}  =  [L1{Aq}  ,  {Au'}  =  lL‘]{Aq'}  ,  {Au"}  =  [L"]{Aq"}  (4.22) 

where  {Aq},  {Aq'}  and  {Aq"}  are  the  nodal  displacement  increment  vectors,  and  [L®],  [L®]  and 
[L®’"]  are  the  corresponding  interpolation  matrices.  Details  of  the  solution  process  are  provided  in 
[22,  23]. 


4.4.3  Coupling  asymptotic  homogenization  with  VCFEM 

The  incremental  energy  functional  for  each  Voronoi  cell  element  in  equation  (4.21)  is  modified  for 
the  asymptotic  homogenization  process  as: 


nf  =  -  /  lst.,,Aa^,-AahdY  -  [  :  A<r^  dY  + 

(o-®  -1-  Aa-y  .  n®  •  (u^  -1-  Au^)  ddY  -  -I-  Ao-®’”  -  <r®®  -  A<r®®)  •  n®  •  (u^'  -|-  Au^')  dY 

-  (<T®®  -f  A<t®®)  •  n®’-  •  (u^"  -h  Au^")  dY-\-  J^(e  +  Ae)A<r^dY  (4.23) 


where  is  an  instantaneous  elastic-plastic  compliance  tensor.  The  last  term  in  equation  (4.23) 
incorporates  the  effect  of  macroscopic  strains  e  in  the  microstructure.  The  stationary  conation 
of  with  respect  to  stress  increment  Aajj  yields  as  Euler’s  equations,  the  incremental  form  of 
kinematic  relations 


e®j  -I-  Aejj  =  e,j  -b  Ae,j  -f 


d(u}  -b  Au}) 
dyj 


(4.24) 


where  Auj  is  the  microscopic  displacement  in  equation  (4.17).  Stationarity  of  the  total  energy 
functional  11  =  11^  with  respect  to  displacement  increments  AnJ,  Auj'  and  Auj",  result 

in  the  inter-element,  interface  and  crack  boundary  traction  reciprocity  conditions  respectively. 
The  Voronoi  cell  finite  element  module  is  incorporated  in  a  macroscopic  analysis  module  with  the 
interface  being  created  by  the  homogenization  procedure  in  [26,  27].  A  displacement  based  finite 
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element  code  with  plane  strain  QUAD4  elements  with  one-point  reduced  integration  and  hourglass 
control  is  used  for  macroscopic  analysis.  Material  constitutive  relations  at  each  integration  point 
of  elements  are  obtained  from  homogenization  of  microscopic  VCFEM  results.  Microscopic  stresses 
tr^are  averaged  to  yield  macroscopic  stresses  S,j.  The  microscopic  VCFEM  is  also  invoked  to 
evaluate  the  homogenized  elastic-plastic  tangent  modulus  by  applying  unit  components  of 
macroscopic  strain. 


4.5  Computational  Subdomain  Level-0  in  the  Hierarchical  Model 

Level-0  corresponds  to  macroscopic  regions  (fi;o  C  ^mat(p))  figure  4.2a,  where  deformation  vari¬ 
ables  like  stresses  and  strains  are  relatively  uniform  in  their  macroscopic  behavior.  Scale  effects  are 
negligible  and  local  constitutive  relations  may  be  derived  from  postulates  of  zero  RVE  volume  in  the 
limit  and  periodicity.  It  is  assumed  that  macroscopic  analyses  with  homogenized  constitutive  rela¬ 
tions  are  sufficient  for  these  regions.  Anisotropic  constitutive  relations  with  varying  parameters  are 
developed  for  continuum  analysis  of  heterogeneous  microstructures  with  elastic-plastic  constituent 
phases.  To  account  for  details  of  microstructural  morphology,  the  constitutive  model  is  based  on 
two  scale  analysis  using  the  asymptotic  homogenization  method  and  microstructural  analysis  byn 
VCFEM.  A  continuum  constitutive  model  can  greatly  enhance  computational  efficiency  over  two- 
scale  analysis. 


4.5.1  An  elastic-plastic  constitutive  model 

Various  continuum  constitutive  models  have  been  proposed,  based  on  unit  ceU  analyses  of  com¬ 
posite  and  porous  microstructures.  One  parameter  plastic  potential  functions  with  assumptions 
of  anisotropy  have  been  introduced  in  [32,  33]  for  composite  materials,  where  the  parameter  is 
determined  by  least  squares  fitting  of  unit  ceU  characteristic  responses.  Bao  et.  al.  [34]  have 
used  the  same  hardening  exponent  for  the  composite  as  for  the  matrix  material.  A  widely  used 
continuum  constitutive  model  for  porous  materials  is  that  of  Gurson  [35],  which  has  been  modified 
by  Tvergaard  [36]  with  unit  ceU  analysis  to  incorporate  the  effects  of  void  growth  and  coalescence. 
Besides  the  Umitations  in  representing  actual  microstructural  heterogeneities,  a  number  of  these 
constitutive  models  do  not  adequately  accommodate  variations  in  constitutive  parameters  with 
evolving  deformation  and  do  not  account  for  post-yield  anisotropy.  Terada  and  Kikuchi  [37]  have 
tried  to  overcome  this  by  using  the  asymptotic  homogenization  to  develop  an  extensive  numerical 
response  database  in  the  strain  space.  Instantaneous  overaU  composite  properties  are  determined 
from  discrete  values  of  homogenized  stress-strain  values  at  points  of  this  database.  This  approach, 
however  leads  to  huge  database  to  cover  aU  possible  deformation  paths  and  requires  solving  an 
inordinately  large  number  of  RVE  boundary  value  problems.  Fish  et.  al.  [10]  have  used  the  idea 
of  transformation  strain  fields,  introduced  by  Dvorak  et.al  [38],  to  develop  a  two  point  averaging 
scheme  based  on  the  mathematical  homogenization  theory  with  piecewise  constant  transformation 
fields.  However,  approximating  the  eigenstrains  with  low  order  polynomial  functions  may  not  be 
able  to  fully  account  for  large  gradients  in  stresses  and  strains  between  phases. 

Motivated  by  two  considerations,  a  piecewise  continuous  elastic-plastic  constitutive  model  with 
an  anisotropic  yield  function  is  developed.  The  first  is  an  accuracy  consideration,  in  that  it  should 
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account  for  the  microstructural  morphology,  e.g.  spatial  distributions,  shapes,  sizes  and  properties 
of  the  individual  phases,  phase  interactions,  as  well  as  the  evolving  stress  and  strain  fields.  This  can 
be  achieved  if  the  model  is  developed  from  detailed  finite  element  analyses  of  the  RVE  (e.g.  VCFEM 
analysis),  subjected  to  a  wide  variety  of  loading  conditions.  The  second  is  an  efficiency  consider¬ 
ation,  since  the  creation  of  a  prohibitively  large  numerical  database  with  a  very  large  number  of 
numerical  experiments  is  of  no  consequence.  The  efficient  development  of  a  constitutive  model, 
accounting  for  underlying  evolution  of  state  variables,  is  accomplished  by  generating  piecewise  con¬ 
tinuous  model  parameters  from  data  in  a  discretized  strain  space  (see  figure  4.4).  Numerical  data 
points  in  the  strain  space  are  systematically  created  through  a  sequence  of  computational  RVE 
analyses  subject  to  an  ordered  set  of  macroscopic  strains  and  strain  paths.  The  strain  space  in 
figure  4.4  is  discretized  into  cubic  elements,  each  containing  32  nodes  or  data  points.  From  the 
computational  RVE  analysis,  constitutive  parameters  like  yield  function  coefficients  and  plastic 
work  are  generated  for  each  nodal  point.  The  constitutive  relation  at  any  point  in  the  strain  space 
are  then  obtained  by  interpolating  nodal  values  using  conventional  shape  functions.  Elastic-plastic 
models  developed  in  the  ensuing  sections  are  for  plane  strain  assumptions. 


Linear  Elasticity 

Orthotropic  homogenized  elastic  material  properties  are  obtained  by  asymptotic  homogenization 
in  conjunction  with  the  VCFE  analysis  of  the  composite  and  porous  microstructures  from  equation 
(4.19)  as  explained  in  [27].  With  plane  strain  assumptions,  three  separate  VCFE  analyses  are 
conducted,  each  corresponding  to  an  independent  component  of  the  macroscopic  strain  {cxx,  ^yy, 
Cxy}.  The  orthotropic  elasticity  tensor  is  stored  for  macroscopic  analysis. 


4.5.2  Elasto-plasticity  with  anisotropic  yield  function 

The  inclusion  phase  in  composites  are  assumed  to  be  linear  elastic,  while  the  matrix  phase  is 
assumed  elastic-plastic  for  both  composite  and  porous  materials.  In  plane  strain  modeling,  an 
assumption  that  the  total  plastic  strain  in  the  out  of  plane  or  ‘third’  direction  is  zero,  is  made. 
The  yield  function  can  then  be  described  in  terms  of  the  macroscopic  in-plane  stress  components 
(Sai,  Syj,  and  S^j,).  The  yield  function  for  porous  and  composite  are  written  using  Hill’s  1948 
anisotropic  yield  function  [39,  33]  in  conjunction  with  the  hydrostatic  stress  dependent  models  of 
the  Tvergaard-Gurson  [35,  36]  as  ; 


$  = 


C(s..-s„)^  +  3s:^ 

y/(w,) 


4-  H  cosh 


( -y/S  (Sjj;  -f  Eyy)\ 

V  2  YfiW,)  ) 


-1  =  0 


(4.25) 


where  exx,eyy,exy  are  the  macroscopic  in- plane  strains,  C{exx,eyy,exy,Wp)  is  a  strain  dependent 
yield  surface  parameter  and  Yf{Wp)  is  the  flow  stress  in  shear.  For  the  composite  materials  the 
dependence  of  pressure  on  yielding  is  deemed  negligible  and  the  hydrostatic  stress  coefficient  H 
is  ignored.  Coefficients  C{exx,€yy,exy,Wp),  Yf(Wp)  and  H  are  determined  from  computational 
experiments  detailed  next.  The  increment  of  plastic  strain  is  obtained  from  the  yield  function  $ 
by  using  the  associated  flow  rule  for  hardening  materials,  i.e. 
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Figure  4.4:  A  Cxx  —  ^yy  —  ^xy  strain  sub-space,  discretized  into  cubic  elements  for  interpolating 
constitutive  model  parameters.  The  nodal  values  of  stresses  and  plastic  work  are  numerically 
generated. 


CoeflRcient  Evaluation 
A.  H  and  Yf(Wp) 


Computational  exercises  indicate  that  the  variation  of  H  with  increasing  hydrostatic  loading  is 
not  significant.  It  is  therefore  assumed  to  be  a  constant  for  all  load  histories.  This  assumption 
is  consistent  with  the  Tverggard-Gurson  models  [35,  36],  where  H  is  determined  in  terms  of  the 
initial  void  volume  fraction.  The  constitutive  parameters  H  and  Yf{Wp)  are  evaluated  in  a  coupled 
manner  by  solving  the  microstructural  EVE  boundary  value  problems  with  two  distinct  loading 
conditions  viz.  (i)  biaxial  tension  loading  (En  =  Syy  =  I^hyd,  ^xy  =  0)  and  (ii)  pure  shear  loading 
(Ejij,  =  Ej/i,  Ea;a;  =  Eyy  =  O).  For  load  condition  (i),  equation  (4.25)  becomes: 

H^hyd.'^hyd,0,Wp)  =  Hcosh  -1  =  0  (4.26) 

and  for  load  condition  (ii),  it  becomes: 

«(0. 0, 2.J,  W,)  =  ^0-^  +  -  1  =  0  or  (4.27) 

The  values  of  H  and  Yf{Wp)  are  determined  iteratively  from  equation  (4.27)  and  further  validated 
against  equation  (4.26).  The  steps  are  as  follows. 

1.  Solve  a  macro- micro  boundary  value  problem  with  EVE  homogenization,  with  incremental 

pure  sheax  loading.  Obtain  macroscopic  plastic  work  by  averaging  the  the  microstructural 
plastic  work  (  Wp  —  Jq  )  and  plot  the  macroscopic  shear  stress  as  a  function  of 

plastic  work  Wp. 

2.  Assume  a  starting  value  for  H  (e.g.  3*/o  as  in  [35,  36])  and  evaluate  Yj{Wp)  from  equation 
(4.27). 

3.  Solve  a  pure  macroscopic  boundary  value  problem  with  incremental  biaxial  loading,  using  the 
homogenized  elastic-plastic  constitutive  relation  and  associated  flow  rule  with  yield  function 
(4.25). 

4.  Plot  the  Cxx  —  ^hyd  and  Cyy  -  J^hyd  curves  for  the  entire  history  of  biaxial  loading. 

5.  Solve  a  macro- micro  boundary  value  problem  with  EVE  homogenization  with  the  same  in¬ 
cremental  biaxial  loading  as  in  the  previous  step.  Plot  the  Y,hyd  —  ^xx  and  T,hyd  —  ^yy  curves 
for  the  entire  history  of  loading. 

6.  Compare  the  results  of  steps  4  and  5.  If  the  two  curves  from  both  methods  are  within  a  preset 
tolerance  everywhere,  then  the  process  is  terminated  and  value  of  H  in  step  2  is  accepted. 
Otherwise,  the  entire  sequence  is  repeated  with  a  dififerent  value  of  H. 

B.  C{eii,Wp) 

The  coefficient  C  in  equation  (4.25)  is  found  to  vary  considerably  with  evolution  of  plastic  de¬ 
formation  and  examples  of  its  variation  with  straining  and  plastic  work  are  shown  in  figures  4.5. 
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(a)  (b) 

Figure  4.5:  Variation  of  the  yield  function  coefficient  C{exx'>€^yy’i^xy<,Wp)^  (a)  as  a  function  of  6  or 
normal  strains  for  different  microstructures,  (b)  as  a  function  of  plastic  work  Wp, 


While  it  is  assumed  to  be  a  function  of  the  total  strain  and  plastic  work,  its  dependence  on  load 
history  is  assumed  to  be  negligible.  In  the  discretized  strain  space  of  figure  4.4,  the  value  of  a  piece- 
wise  continuous  C  at  any  point  may  then  be  obtained  by  interpolation  from  nodal  values  according 
to 

32 

9  ^yy^  ^xy^  Wp)  =  E  C'a(Wp)iV„(  ^xx^  ^yy  5  ^xy  )  (4.28) 

a=l 

where  Ca  are  the  nodal  values  and  No,  are  shape  functions  for  a  32-noded  brick  element. 


Generation  of  Discretized  Strain  Space  and  Nodal  Parameters 

The  nodal  values  of  macroscopic  stresses  (Sx®?  ^yy,  Sa;y)  and  the  corresponding  plastic  work  Wp  are 
first  evaluated  at  each  nodal  point  in  a  subspace  of  the  e®®  —  eyy  —  e^y  space  by  solving  incremental 
macro-microscopic  boundary  value  problems  with  VCFEM  and  asymptotic  homogenization.  In 
this  process,  macroscopic  strain  increments  are  applied  to  the  RVE  subjected  to  periodic  boundary 
conditions  [26,  27].  Strain  increments  are  applied  along  a  radial  line  in  the  strain  space,  such  that 
a  constant  ratio  between  strain  components  is  maintained,  i.e.  Ae®®  :  Aej,j,  :  Ae®y  =  1  :  tanO  : 
(1  -f-  tan^6)tan<f>,  where  6  and  (p  are  the  angular  coordinates  in  the  strain  space  of  figure  4.4.  The 
fiow  stress  Y}{Wp)  at  each  node  in  figure  4.4  can  be  obtained  from  the  shear  stress-plastic  work  plot 
and  equation  (4.27).  From  the  values  of  macroscopic  variables  (E®®,  Eyy,S®y,y/(Wp))  at  a  node 
corresponding  to  the  end  of  an  increment,  the  coefficient  C'(e®®,eyy,e®y,  Wp)  is  calculated  using 
equation  (4.25). 
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From  the  symmetry  conditions,  only  a  quarter  of  the  e^x  —  ^yy  —  ^xy  strain  space  is  considered 
for  loading  such  that  Q°  <6  <  180°  in  the  Cxx  —  tyy  plane  and  0°  <  ^  <  90°  outside  of  this  plane. 
This  is  indicated  in  figure  4.4.  This  chosen  subspace  of  the  strain  space  is  divided  into  16  cubic 
brick  element  with  32  nodal  points  each.  The  location  of  elements  in  the  strain  space  are  selected 
to  optimally  account  for  the  variations  in  C.  These  variations  in  C  with  the  coordinate  angle 
0  (location  in  the  Cxx  —  ^yy  plane)  and  the  plastic  work  Wp  for  the  different  microstructures  are 
plotted  in  figure  4.5.  The  parabolic  form  of  C  in  figure  4.5a  is  consistent  with  the  quadratic  term 
(^xx  ^yy  in  the  yield  function.  The  minimum  values  occur  near  6  =  135°  corresponding  to  a 
pure  deviatoric  state.  The  coeflRcient  subsequently  increases  to  account  for  the  increase  in  plastic 
work  in  the  yield  function  $.  In  figure  4.5b,  the  coefficient  C  as  a  function  of  the  plastic  work, 
which  corresponds  to  the  radial  direction  in  the  strain  space,  is  plotted.  The  value  of  C  stabilizes 
beyond  a  value  of  the  plastic  work,  which  is  used  as  the  outer  boundary  of  the  strain  space  envelope 
in  figure  4.4. 


4.5.3  Numerical  implementation  of  the  constitutive  model 

The  elastic-plastic  constitutive  model  for  composite  and  porous  materials  is  derived  from  the 
anisotropic  yield  function  (4.25)  with  associated  flow  rule  and  isotropic  hardening.  In  an  incre¬ 
mental  form,  the  stress  increments  AS,j  are  related  to  elastic  increments  of  strains  (Ae^/  —  Ae^^) 
admitting  additive  decomposition,  as 

ASij  =  £§fc,(Aefc/  -  Ae",)  (4.29) 

where  is  the  homogenized  elasticity  tensor.  Using  associated  flow  rule,  components  of  the 
plastic  strain  increment  are  obtained  as: 


Aef.  =  AA 


d'Ziij 


(4.30) 


Elimination  of  the  flow  parameter  AA  from  the  above  equations  results  in  the  two  equations 


These  equations  are  solved  using  the  backward  Euler  integration  method,  with  gradients  evaluated 
at  the  end  of  the  increment.  With  known  increments  of  strain,  the  resulting  set  of  equations  (4.31) 
together  with  the  yield  function  (4.25)  are  solved  iteratively  by  using  the  Newton- Raphson  method. 
The  stress  increments  are  obtained  by  the  following  steps. 

1.  Initialize  values  of  AS^x,  ASj^j,  and  AE^y. 

2.  Calculate  the  gradient  (^■)  of  the  yield  function  and  solve  for  the  increments  of  plastic 
strain  AeJ^,,  Ae^^^  and  Aegj,  from  equations  (4.29)  and  (4.31).  Update  the  stresses  and 
plastic  work  using  the  relation  AWp  =  S^xAeP^,  -1-  Syy  Ae^j^  +  S^yAe^j,. 

3.  If  $  <  toll  and  correction  to  plastic  strain  increment  Se^j  <  to/2,  where  toll  s-^d  to/2  are 
prescribed  tolerances,  then  stop.  Otherwise  go  to  step  2. 

The  parameter  C  is  then  obtained  from  the  interpolation  equation  (4.28). 
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4.5.4  Numerical  examples  with  the  continuum  constitutive  model 

The  elastic- plastic  constitutive  model  is  validated  by  comparing  the  results  of  macroscopic  numerical 
simulations  with  those  generated  by  macro-micro  scale  analysis  using  asymptotic  homogenization. 
Examples  are  conducted  for  both  composite  and  porous  materials  with  different  microstructure 
morphologies,  viz.  different  shapes,  sizes  and  spatial  distributions  of  the  heterogeneities. 

Analysis  of  Composite  Microstructures 

Six  microstructural  RVE’s  of  20%  Vj  Alumina-aluminum  composite  with  Alumina  fiber  in  alu- 
nainum  matrix,  as  shown  in  figure  4.6,  are  analyzed.  The  RVEs  are  classified  as: 

(a) :  Square  edge  pattern  with  a  circular  inclusion  (Cl) 

(b) :  Square  edge  pattern  with  an  elliptical  inclusion  (aspect  ratio  f  =  3)  (C2) 

(c) :  Random  pattern  with  25  identical  circular  inclusions  (C3) 

(d) :  Horizontally  aligned  random  pattern  with  25  identical  elliptical  inclusions  (C4) 

(e) :  Randomly  oriented  random  pattern  with  25  identical  elliptical  inclusions  (C5) 

(f) :  Random  pattern  with  17  random  shape  and  size  inclusions  (C6) 

The  material  properties  for  the  elastic  Alumina  fiber  are: 

Young’s  Modulus  {Ec)=  344.5GPa,  Poisson  Ratio  (^c)=  0.26; 
and  for  the  elastic-plastic  Aluminum  matrix  are: 

Young’s  Modulus  (Em)=  68.3  GPa,  Poisson  Ratio  (i/m)=0.30.  Initial  Yield  Stress  (lb):  55  MPa, 
and  Post  Yield  hardening  law:  aggv  =  lb  +  2.08€Pj^. 

The  RVEs  are  subjected  to  four  different  types  of  loading  viz. 

•  LI:  Pure  shear  loading  with  increments(ASj;j;  =  AEyy  =  0,  =  AS) 

•  L2:  Uniaxial  tension  loading  with  increments  (AS^jj;  =  AS,  AEyy  =  AH^y  =  0) 

•  L3:  Biaxial  tension  loading  with  increments  (ASa,®  =  ASj/j,  =  AS,  AS^j,  =  0) 

•  L4:  Biaxial  tension-compression  loading  with  increments  (AS^®  =  -A^yy  =  AS,  AS®^  =  0) 

The  stress-strain  response  of  the  six  composite  microstructural  RVEs  with  the  four  loading 
conditions  are  conducted  and  the  results  for  simple  tension  (LI)  are  depicted  in  figure  4.7.  The 
results  by  the  constitutive  model  and  two-scale  asymptotic  homogenization  approach  are  generally 
found  to  agree  very  well  for  the  entire  range  of  loading  upto  fairly  high  level  of  straining.  The 
only  discrepancy  is  found  with  the  biaxial  tension  loading  condition  (L3),  for  which  the  devia¬ 
tion  strains  are  shown  in  table  1.  However  the  deviations  occurs  at  high  strains  levels,  for  which 
the  stresses  are  nearly  twice  the  matrix  yield  stress.  It  is  important  to  note  that  the  deviation 
of  continuum  model  response  from  the  two-scale  asymptotic  homogenization  response  can  be  used 
as  a  signal  for  switching  from  the  former  to  the  latter  type  of  analysis  in  a  multiple  scale  simulation. 

As  an  example  of  a  structural  analysis  with  the  two  approaches,  a  square  plate  with  a  square 
hole  is  solved  with  tension  loading  on  two  opposite  faces.  A  quarter  of  the  plate  with  symmetry 
and  loading  conditions  is  shown  in  figure  4.8a.  A  total  traction  of  55  MPa  is  applied  in  10  equal 
increments.  The  material  is  a  20  %  Vj  Alumina- aluminum  composite  with  the  microstructural 
RVE  consisting  of  15  identical  circular  inclusions  dispersed  randomly  in  the  matrix  (figure  4.8b). 
The  material  properties  are  the  same  as  in  the  previous  example.  Comparison  of  the  results  are 
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(Cl, VI) 


(C2,V2) 


(C3,V3) 


(C4,V4) 


(C5,V5) 


(C6,V6) 


Figure  4.6:  Microstructures  with  different  shape,  size,  orientation  and  spatial  distribution,  for  20% 
volume  fraction  composite  (C1,C2,C3,C4,C5,C6)  and  porous  (V1,V2,V3,V4,V5,V6)  materials. 
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(a) 


axial  atrain  in  loading  diiwction 


(b) 


(c) 


(d) 


axial  strain  In  loading  diroction 

(f) 


Figure  4.7:  Comparison  of  the  stress-strain  results  in  the  composite  for  the  uniaxial  loading  test 
by  (i)  macroscopic  analysis  with  the  homogenized  constitutive  model  and  (ii)  two-scale  analysis 
with  asymptotic  homogenization  model  for:  (aj^l,  (b)  C2,  (c)  C3,  (d)  C4,  (e)  C5  and  (f)  C6 
microstructures. 


iririririy 


Figure  4.8:  (a)  Finite  element  model  for  a  quarter  of  the  square  composite  plate  with  square  hole 
and  (b)  the  VCFE  model  of  the  microstructural  composite  RVE. 


made  through  contour  plots  of  the  macroscopic  stress  Hxx  (shown  in  figure  4.9)  and  macroscopic 
plastic  work  Wp  (not  shown).  The  figures  reveal  that  at  all  locations  the  difference  in  the  two 
approaches  is  less  than  1%.  However,  the  computational  efficiency  of  the  macroscopic  analysis 
with  the  continuum  constitutive  model  is  far  superior  than  the  two  scale  analysis.  The  scale  up 
in  efficiency  for  this  problem  is  approximately  75000%,  and  is  therefore  very  desirable  when  only 
macroscopic  results  are  of  interest. 

Analysis  of  Porous  Microstructures 

The  six  microstructural  RVEs  are  analyzed  again  for  porous  materials,  with  voids  replacing  the 
inclusions  in  figure  4.6.  The  material  considered  is  an  aluminum  alloy  with  20%  void  volume 
fraction  and  the  following  elastic-plastic  properties: 

Young’s  Modulus  (Em)=  68.3  GPa,  Poisson  Ratio  0.30,  Initial  Yield  Stress  (ib)=  55  MPa, 

and  Post  Yield  hardening  law:  Ceqv  =  Fb  +  2.08e§^„.  The  same  four  load  histories  (L1,L2,L3,L4) 
are  applied.  An  important  difference  between  the  composite  and  porous  microstructures,  is  that 
in  the  latter  plastic  strain  localization  in  small  regions  is  a  common  occurrence  depending  on  the 
void  morphology  and  the  nature  of  loading.  Such  nonhomogeneous  distribution  of  plastic  strains 
is  a  major  source  of  discrepancies  between  responses  by  the  two  approaches  and  act  as  ‘limiters’ 
for  the  range  of  application  of  the  continuum  model.  Microstructures  with  homogeneous  and 
nonhomogeneous  distributions  of  plastic  strain  are  shown  in  figure  4.10.  For  the  microstructure 
VI  (square  edge  pattern  with  a  circular  void)  the  strain  distribution  is  quite  uniform  in  pure 
shear  loading,  while  for  the  microstructure  V2  (square  edge  with  an  elliptical  void)  there  is  intense 
localization  with  narrow  ligaments.  Consequently  the  continuum  model  ceases  to  be  effective. 

As  in  the  case  with  composites,  the  main  challenge  for  the  homogenized  constitutive  model 
is  encountered  during  simulations  with  biaxial  tension  loading,  i.e.  (ASi,®  =  =  ASj  and 
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Max. 


1.219E-01 


Min. 


1.047E'01 


e.746E-02 


7.027E-02 


5.307E-02 


3.588E-02 


(a)  (b) 

Figure  4.9:  Comparison  of  macroscopic  stress  {'Exx)  contours  in  composite  material  with  (a)  the 
homogenized  constitutive  model  and  (b)  two-scale  analysis  with  asymptotic  homogenization. 

ASjjj,  =  0).  The  first  term  in  the  yield  function  (4.25)  drops  out  for  this  loading  and  the  model 
delivers  the  same  amounts  of  plastic  strains  in  the  x  and  y  directions.  Due  to  the  lack  of  anisotropy 
in  the  hydrostatic  term  in  the  yield  function,  the  continuum  model  is  effective  only  for  those 
microstructures  that  exhibit  near-isotropic  plastic  behavior  for  this  loading.  Microstructures  V2 
(square  edge  with  elliptical  void)  and  V4  (horizontally  aligned  elliptical  voids  with  random  spatial 
distribution)  shows  very  different  strain  responses  for  each  direction  with  biaxial  loading.  Thus  the 
continuum  constitutive  model  is  largely  ineffective  for  these  RVEs.  The  microstructures  V5  (ran¬ 
domly  oriented  identical  elliptical  voids  with  random  spatial  distribution)  and  V6  (random  spatial 
distribution  with  random  shape  and  size)  also  show  significant  plastically  induced  directional  effects 
and  the  constitutive  models  are  therefore  restricted  to  the  elastic  range.  The  list  of  performance 
and  strain  ranges  of  all  the  microstructures  with  the  different  loadings  are  given  in  table  2.  The 
microstructures  VI  and  V3  exhibit  relatively  isotropic  responses  and  yield  satisfactory  agreement 
between  responses  by  the  constitutive  model  and  the  two-scale  asymptotic  homogenization.  Com¬ 
parisons  of  stress-strain  responses  for  biaxially  and  uniaxiaUy  loaded  RVEs  are  made  in  figures  4.11. 
These  show  very  good  agreement.  The  RVEs  V2  and  V4  exhibit  intense  localization  early  on  in 
the  straining,  while  the  RVEs  V5  and  V6  exhibit  marked  anisotropy  with  biaxial  loading.  Plastic 
flow  predictions  for  these  RVE’s  with  the  continuum  model  are  therefore  not  reliable. 

4.6  Computational  Subdomain  Level- 1  in  the  Hierarchical  Model 

The  level-1  regions  (D;i  C  ^mat{p))  the  computational  domain  are  identified  by  locally  high 
gradients  of  macroscopic  variables  e.g.  stresses,  strains  strain  energy  etc..  These  are  generally  in¬ 
dicators  of  imminent  damage  evolution  or  localization  in  the  microstructure.  In  these  regions,  the 
assumptions  of  macroscopic  uniformity  and  statistical  periodicity  of  the  RVE  are  stiU  assumed  to 
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(b) 

Figure  4.10:  Contour  plot  of  the  microscopic  plastic  strain  in  the  voided  microstructures  under 
pure  shear  loading  condition,  for  (a)  VI  (b)  V2  microstructures. 
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(a) 


(b) 


(c) 


(d) 


Figure  4.11:  Comparison  of  the  stress-strain  results  in  the  porous  material  for  (a)  bi-axial  ten¬ 
sion  for  VI  microstructure  (b)  bi-axial  tension  for  V3  microstructure,  (c)  uniaxial  tension  for  VI 
microstructure  and  (d)  uniaxial  tension  for  V3  microstructure. 
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Figure  4.12:  Creating  extended  microstructure;  (a)  a  mesh  of  macroscopic  elements  with  an  under¬ 
lying  microstructure  of  repeated  RVEs,  (b)  the  extended  microstructure  by  Voronoi  tessellation. 


hold.  However,  these  assumptions  may  cease  to  hold  soon,  with  additional  deformation.  Thus,  in 
concurrence  with  macroscopic  simulations,  computations  are  necessary  at  the  microstructural  scale 
to  monitor  the  initiation  and  growth  of  damage  and  instabilities.  The  RVE  response  for  level- 1 
elements  is  explicitly  calculated  using  asymptotic  homogenization  and  VCFEM  analysis.  In  its 
computational  implementation,  a  sequence  of  finite  element  analyses  are  executed  as  follows: 

(a)  Microstructural  analysis  by  VCFEM  with  unit  strains  or  increments  and  periodicity  boundary 
conditions  to  generate  homogenized  material  parameters 

(b)  Macroscopic  analysis  of  structural  components  with  homogenized  coeificients, 

(c)  Microstructural  analysis  by  VCFEM  with  imposed  macroscopic  strains  or  increments  and  period¬ 
icity  boundary  conditions  for  simulating  the  evolution  of  microscopic  stresses,  strains  (o-',  e')  etc.  . 
The  computational  requirements  of  this  level  are  considerably  higher  than  that  for  level  0,  since 
at  each  integration  point  in  the  macroscopic  computational  mesh,  a  complete  microstructural  RVE 
problem  has  to  be  solved  twice.  It  is  also  critical  that  this  level  be  discontinued  once  the  microstruc¬ 
tural  damage  or  instability  evolved  beyond  pre- determined  threshold  values. 


4.7  Computational  Subdomain  Level-2  in  the  Hierarchical  Model 

Level-2  regions  (figure  4.2b)  are  classified  as  those  with  severe  microstructural  nonuniformities  in 
the  form  of  evolving  damage.  This  results  in  loss  of  statistical  periodicity  of  the  assumed  RVE  and 
these  regions  may  be  identified  with  ^MAT{np)  in  equation  (4.5).  In  the  computational  model,  the 
level-2  elements  (fi;2  C  ^MAT{np))  materialize  from  the  microstructure  of  level-1  elements.  It  is 
assumed  that  prior  to  this  transition,  the  level-1  elements  have  been  refined  to  reach  sufficiently  high 
spatial  resolution.  In  the  macro-micro  model  of  level- 1  switches  to  a  completely  microscopic 
model. 
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The  method  of  generating  the  level-2  microstructure  in  each  element  is  iUustrated  in  figure 
(4.12).  An  extended  microstructure  is  first  created  by  repeating  RVEs  in  succession,  to  cover  the 
entire  region  of  the  macroscopic  level-1  elements  in  transition  to  level-2.  In  equation  (4.7),  a  local 
non-periodic  region  is  created  as  : 

^MAT(np)  ^  (4  32) 

where  Y'^  corresponds  to  the  RVEs  in  a  periodic  domain,  adjusted  for  microscopic  evolution.  This 
is  then  overlayed  by  the  macroscopic  elements  to  accurately  encase  the  level-2  region  with  clearly 
delineated  boundaries.  Each  level-2  element  now  contains  a  heterogeneous  material  distribution 
(Y|)  that  is  defined  as  the  intersection  of  the  non-periodic  material  region  and  the 

element  domain  llfj 

=  Yf  n  Of2  (4.33) 

This  region  is  subsequently  tessellated  into  a  mesh  of  Voronoi  ceU  elements  (figure  4.12b).  Trac¬ 
tion  continuity  between  level-2  microstructural  region  and  neighboring  level-O/level-1  elements  is 
incorporated  through  a  layer  of  transition  elements. 

4.7.1  Damage  criterion  for  particle  cracking  in  level-2 

The  level-2  VCFEM  modeling  consist  of  brittle  reinforcing  particles  and  a  ductile  matrix  material. 
For  the  brittle  particulate  materials,  microstructural  damage  initiation  is  assumed  to  be  governed 
by  a  maximum  principal  stress  based  criterion.  In  this  criterion,  a  crack  is  initiated  when  the 
maximum  principal  stress  in  tension  exceeds  a  critical  fracture  stress  tr"'  at  a  point.  In  the  compu¬ 
tational  procedure,  complete  particle  splitting  is  assumed  to  occur  in  the  form  of  an  elliptical  void, 
normal  to  the  principal  direction,  as  soon  as  the  principal  tensile  stress  reaches  In  the  case 
of  particle  splitting,  the  crack  tip  extends  nominally  into  the  matrix.  In  the  incremental  compu¬ 
tational  procedure,  more  than  one  point  may  exceed  the  critical  <7‘^  value  during  increment.  The 
location  of  a  single  crack  is  determined  by  a  weighted  averaging  method  as: 

^damage  =  <T’f{x,y)  ’  Vdamage  =  iT^x,y)  V)  ^  (4-34) 

where  crj{x,y)  corresponds  to  all  values  of  maximum  tensile  principal  stress  larger  than  (7‘^  in  the 
particle.  The  crack  is  oriented  at  right  angles  to  the  principal  stress  directions  at  {x damage,  Vdamage) 
and  extends  to  the  interface  on  both  sides. 


4.8  Coupling  the  Levels  in  the  Hierarchical  FEM 

While  level-0  elements  (E/o  €  fl/o)  and  level-1  elements  (E/i  €  fl;i)  are  coupled  naturally  through 
the  familiar  assembly  process,  the  interfacing  of  level-2  (E/2  €  11/2 )  elements  with  either  of  the  first 
two  requires  more  attention.  The  mismatch  in  the  number  of  boundary  nodes  in  these  elements 
necessitate  the  introduction  of  transition  elements  (Et,  G  ^12)1  acting  as  buffer  zones  as  shown  in 
figure  4.13.  Both  E/2  and  Etr  elements  employ  VCFEM  for  setting  up  the  element  stiffness  matrices 
and  load  vectors.  The  entire  computational  domain  is  thus  comprised  of 

=  {ft/o  U  fi/i  U  fl/2  :  fi/o  =  u£?iE/o;  fi/i  =  u£\E/i;  fi/2  =  u£\E/2  U 
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Figure  4.13:  Interfaces  between  the  levelO/1  Eio  and  Ei\  elements,  transition  Etr  elements  and  the 
level-2  Eli  elements. 


for  which  the  nodes  are  differentiated  as  (see  figure  4.13): 

(i)  or  internal  level-2  Voronoi  cell  element  nodes  that  are  not  on  any  interface  or  boundary, 

(ii)  (nd'2  '  or  Voronoi  cell  element  nodes  on  the  E12-E12  interface, 

(iii)  (nd'2/<’-  )  or  Voronoi  ceU  element  nodes  on  En-Etr  interface, 

(iv)  (nd^oi/tr 

)  or  nodes  on  Eioji-Etr  interface  that  belong  to  Eio  and  En  elements, 

(v)  (nd**^)  or  Voronoi  ceU  element  nodes  on  Eioji-Efr  interface  that  do  not  belong  to  Eio  and  En 
and  need  to  be  staticaUy  condensed.  In  an  incremental  analysis  for  elasto-plasticity,  the  principle 
of  virtual  work  for  the  computational  domain  at  the  end  of  the  n  —  th  increment  may  be  written 
as  a  scalar  valued  function  G  in  terms  of  variables  at  different  levels  as: 


(?”+^(Au,«u)  = 


+ 

+ 


[  ^ij{n-  +  An)^da-  f  JiSuidQ 

f  -t- Au)^dQ  -  /  fiSuidO. 

f  +  Au)^^dfi  -  /  JibuidO. 

I  tiSuidT  -  f  tiSuidT  -  f  USuidT 

./r”+^  Jr”^^  Jr”^^ 


(4.35) 


In  an  iterative  solution  process  with  the  Newton- Raphson  method,  a  consistent  linearization  by 
taking  the  directional  derivative  of  along  incremental  displacement  vector  Au,  yields  the 

tangent  stiffness  matrix.  For  the  i  —  th  iteration,  the  assembled  equations  have  the  foUowing 


58 


structure. 


■  J^10/1,10/1  J^lO/1,12 

J^12,10/1  J^12,12 

Here  corresponds  to  displacements  at  nodes  (nd^°/^)  that  belong  to  elements  Eiq  and  En 

in  the  computational  region  U  ft/j  as  shown  in  figure  4.13.  It  should  be  noted  that  they  also 
include  nodes  at  the  interface  of  elements  Eio/En  and  elements  Etr-  The  displacements  A on 
the  other  hand  corresponds  to  nodes  on  the  interfaces  of  elements  Etr  and  elements  E12  or 

to  nodes  (nd*^/^)  on  the  interfaces  between  two  E12  elements.  Contributions  to  the  stiffness  matrix 
\K\  and  load  vector  {F}  for  elements  in  U  fl/j  may  be  obtained  as 


AU^o/i 

^U‘2B 


AF'2 


(4.36) 


JQ 


dN^  dT.ii  dN9 


dxi  deki  dxi 

*0  *1 


dO, 


/  fiN^dn+  f  tiN?dT-  [ 

^  “  yrr+'urr+'  h 


rn+iur"+i 
^0  *1 


•  dN^ 

,  Eij{un  +  Au')^^dQ 

/nr+'ufir+'  dxj 


%  '"“ll 


a,/?  correspond  to  the  node  numbers  and  Na  are  the  shape  functions  of  macroscopic  elements. 
For  the  elements  Etr  and  Ei2  belonging  to  contributions  to  the  stiffness  matrix  and  load 
vectors  are  obtained  by  VCFEM  calculations,  together  with  static  condensation.  The  transition 
element  facilitates  continuous  variation  from  microscopic  to  macroscopic  elements  without  jumps 
or  discontinuities.  On  the  Eio/i-Etr  interfaces,  this  is  accomplished  by  constraining  displacements 
at  nodes  nd®^  to  conform  to  displacement  interpolation  of  the  adjacent  Eio  or  En  elements.  The 
constraints  at  the  nodes  nd®^  are  applied  in  terms  of  displacements  at  as 

{At/*"}  =  [Q]  {At/'®^/*’'}  (4.37) 


where  [Q]  is  the  constraining  matrix.  For  bilinear  QUAD4  level-0/1  elements,  each  row  of  [^] 
consists  of  the  inverse  of  the  distance  of  the  constrained  node  to  the  corner  nodes.  The  interfaces 
with  the  Ei2  elements,  i.e.  the  Ej2-Etr  interfaces  are  treated  in  the  same  way  as  Ei2-Ei2  interfaces. 


The  displacements  in  equation  (4.36)  correspond  to  those  at  the  boundary  nodes  (nd^^/^, 

nd^2/fr,  and  nd*")  of  level-2  elements  that  contain  the  microstructural  VCFE  model.  To 

account  for  the  contribution  of  the  internal  nodes  (nd^^^),  it  is  therefore  necessary  to  use  static 
condensation  and  recovery  process  for  representing  the  VCFEM  displacement  solutions  at  the 
internal  nodes  in  terms  of  the  displacements  at  boundary  nodes  AU^^^,  where 


{  At/'2B  J  ^ 


AJ7/2/2  ' 

AU^Vir 
A[/-Wl/<r  ' 

At/*" 


7  0  0  0 
0/00 
o  0  7  0 
0  0  0  Q 


'  At/^2/2  'I 

^  AI/'2/‘^  I 

^  At/'°1/*’'  J 


(4.38) 


where  [I]  and  [Q]  are  the  identity  and  constraint  matrices  respectively.  In  a  VCFEM  solution 
process  the  stiffness  matrix  and  the  load  vector  may  be  partitioned  accordingly  as 

■  J^12B,12B  J^12B,12I 

J^12I,12B  J^12I,12I 


f  At/'25  \  _f 
1  Al/'2/  /  -  i 


AJ7’/2B 

AF'2/ 


(4.39) 


59 


Static  condensation  of  the  internal  degrees  of  freedom  yields 

_  Jjj'/2B,;2/j  jjj'i2/,/2/j  jj^'/2/,;2Bj  j  |^[//2BJ  _ 

(4.40) 

These  stiffness  matrices  and  load  vectors  are  then  used  in  the  assembly  process  of  equation  (4.36). 


4.9  Adaptation  Criteria  for  Various  Levels 

It  is  evident  that  appropriate  criteria  need  to  be  established  for  transitioning  the  computational 
subdomains  from  one  level  to  another.  In  addition  to  level  transitions,  element  refinement  by  h- 
adaptation  is  also  executed  in  level-0  and  level-1  regions  for  two  reasons.  The  first  is  to  identify  and 
reduce  any  suitably  chosen  ‘error  measure’  in  the  computational  model.  Secondly,  the  h-adaptation 
enables  the  computational  model  to  ‘zoom  in’  on  regions  of  evolving  localization  due  to  microscopic 
non-homogeneity.  It  reduces  the  disparity  in  the  macro-  and  micro-  scale  elements  by  successive 
refinement  of  macroscopic  elements  in  the  critical  regions  as  shown  in  figures  (4.2a,4.16).  The  cri¬ 
teria  used  for  element  refinement  and  level  transition  are  delineated  below. 

•  Level  0/ 1  h-adaptation  may  be  executed  based  on  an  error  measure  may  be  defined  as 

Ee  =  |l/(u  -  Uh)||  (4.41) 

where  u  is  the  true  solution,  is  its  finite  element  approximation,  /  is  any  appropriate  function  of 
deformation  e.g.  energy,  stresses,  strains  etc.,  and  ||  •  ||  denotes  a  norm  of  the  function.  Substantial 
work  on  h-method  of  mesh  refinement  e.g.  [40,  41,  42]  have  proposed  various  forms  of  /  and  norms 
to  optimize  the  computational  effort  and  improve  reliability.  In  this  work,  elements  are  subdivided 
based  on  an  element  level  parameter  defined  as 

/  IP  ^  I|/('^h)ll  fA  /ION 

^  ^  ^ 

where  (Ee)cr  corresponds  to  a  threshold  element  parameter,  above  which  the  element  needs  to  be 
refined  and  qi  is  a  prescribed  quality  index.  ||/(uii)||  =  ||/(uh)l|e  is  the  sum  of  element  norms 

of  any  function  for  the  entire  computational  domain  and  NE  is  the  number  of  elements.  In  this 
work  ||/(uii)||e  is  taken  to  be  the  maximum  difference  in  plastic  work  with  neighboring  elements. 
Accordingly  a  value  of  91=1.5  is  found  to  be  satisfactory  for  this  study. 

•  Level-0  to  level-1  transition  of  Eio  elements  takes  place  in  accordance  with  criteria  based 
on  macroscopic  variables  e.g.  (17,  e)  in  the  continuum  model  that  are  related  to  critical  micro¬ 
scopic  variables.  Additionally,  this  transition  is  activated  when  the  continuum  level  model  in  section 
4.5.1  starts  to  digresses  considerably  from  the  results  of  two-scale  homogenization.  Different  criteria 
are  used  for  composite  and  porous  materials. 

A.  For  composite  microstructure  with  inclusions: 

(E/)  >  S"’- =  r  *  <r"  and/or  ).  >  (4.43) 

Here  S/  is  the  maximum  principal  stress  in  Eio,  S‘^’’  is  a  critical  stress  that  is  a  pre-determined 
fraction  r  of  the  critical  fracture  stress  {(t")  for  microstructural  inclusions.  The  fraction  r  is 
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selected  to  be  |  in  this  study.  The  second  condition  is  established  since  the  earliest  digression  from 
the  homogenization  results  is  observed  for  biaxial  loading,  and  established  from  the  results  of 
section  4.5.4. 

B.  For  porous  microstructure  with  voids:  Strain  based  criteria  are  deemed  more  important  in  the 
case  of  damage  in  porous  materials  and  hence  transition  is  activated  when 

(^)  >  (^r  =  r  *  and/or  (4.44) 

where  e^(=  is  the  macroscopic  effective  plastic  strain  in  Eiq  and  is  a  pre-determined 

fraction  r  of  a  suitably  assessed  high  microstructural  plastic  strain  Again  is  established 

from  the  limiting  values  of  biaxial  strain  in  section  4.5.4. 


•  Level-l  to  level-2  transition  of  Ew  elements  occur  when  a  sufficiently  high  spatial  resolu¬ 
tion  is  attained  by  h-refinement  and  when  the  RVE  is  deemed  to  be  on  the  verge  of  significant 
damage  evolution.  The  adaptive  criteria,  which  monitor  the  transition  from  elements  En  to  E12 
are  prescribed  in  terms  of  topological  evolution  of  microscopic  damage  as: 

A.  For  composite  microstructure  with  inclusions: 


7?.  ~  ^ 

J^damage  — 

B.  For  porous  microstructure  with  voids: 

damage 


_adr 

•^damaoR  —  .  ^  ^ 


AR  - 


(4.45) 


(4.46) 


where  Rdamage  Is  the  ratio  of  the  number  of  damaged  inclusions  (NPD)  to  the  total  number  of 
inclusions  (NP)  for  composites.  For  porous  microstructures,  it  is  the  ratio  of  the  damaged  area 
corresponding  to  regions  which  have  high  plastic  strains  to  the  total  RVE  area.  Rcr  is  a  prescribed 
critical  ratio  and  varies  from  problem  to  problem. 


4.10  Numerical  Examples  with  the  Adaptive  Multilevel  Model 

Numerical  examples  are  solved  to  understand  the  effect  of  structural  geometry  and  microstructural 
morphologies  on  the  macroscopic  and  microscopic  responses.  In  all  examples,  the  inclusions  are 
assumed  to  be  linear  elastic  which  can  crack  by  the  principal  stress  criterion  and  the  matrix  is  elastic- 
plastic.  In  the  first  example,  a  RVE  consisting  of  a  single  circular  inclusion  is  modeled  by  level-1 
and  level-2  elements  under  applied  uniaxial  tension  loading.  In  the  first  case  the  En  macroscopic 
element  is  coupled  with  the  microscopically  periodic  RVE  by  asymptotic  homogenization,  while  in 
the  second  case,  the  tension  load  is  directly  applied  on  one  edge  of  the  Voronoi  cell  element  Ei2  model 
of  the  RVE.  The  loading  is  continued  beyond  the  level  of  inclusion  cracking.  This  is  represented 
by  the  loss  of  stress  carrying  capacity  in  the  averaged  stress-stain  plot  of  figure  4.14a.  Though  the 
response  is  similar  in  the  initial  elastic  phase,  it  diverges  with  the  onset  of  plastic  fiow.  The  damage 
occurs  earlier  in  the  level- 1  element  due  to  the  additional  constraint  of  periodicity.  The  contour 
plots  of  inclusion  maximum  principal  stress,  normalized  by  the  critical  stress  cr"",  are  compared 
at  a  macroscopic  strain  2.1%  in  figure  4.14.  This  is  just  before  cracking  by  the  level-1  model.  A 


61 


considerably  lower  stress  level  is  seen  for  the  level-2  model.  Such  over-estimation  of  stresses  with 
periodicity  boundary  condition  near  free  boundaries  necessitate  the  use  of  the  proposed  multi-level 
models. 


4.10.1  Effects  of  inhomogeneity  size  (Scale  effect) 

Neglecting  scale  effects,  that  reflect  the  actual  size  of  microstructural  constituents,  is  a  charac¬ 
teristic  of  homogenization  with  assumptions  of  statistical  periodicity  of  a  vanishingly  small  RVE. 
While  elements  Eio  and  En  conform  to  this  restriction,  the  level-2  elements  E12  in  this  work  depict 
the  actual  size  of  the  microstructure  through  the  multi-level  coupling  and  hence  scale  effects  prevail. 

In  this  example  the  effect  of  particle  or  void  size  on  the  evolution  of  damage  is  investigated.  The 
two  different  microstructural  RVE’s  considered  in  figure  4.15b  (i  and  ii)  have  identical  distribution 
(square  edge)  and  same  particle  or  void  area  fraction  oiV  —  f  =  20%.  But  the  RVE  sizes  are 
different  in  that,  the  size  of  the  smaller  RVE(i)  is  0.5mm  x  0.5mm  while  that  of  the  larger  RVE(ii) 
is  1mm  X  1mm.  The  macroscopic  structure  is  a  square  plate  with  a  square  hole,  for  which  the 
initial  level-0  mesh  with  dimensions  is  shown  in  figure  4.15a.  Only  a  quarter  of  the  plate  is  modeled 
from  symmetry  considerations.  The  smallest  size  of  macroscopic  element  allowed  in  this  analysis 
by  h  —  adaptation  is  set  to  1mm.  Thus  the  Ei2  elements  contain  only  one  element  in  the  VCFEM 
model  for  the  larger  RVE(ii)  but  four  elements  in  the  VCFEM  model  for  the  smaller  RVE(i).  The 
material  properties  for  both  composite  and  porous  materials  are  as  follows. 

Aluminum  matrix  (Elastic-Plastic);  Young’s  Modulus  (Em)=68.3  GPa,  Poisson  Ratio  (i/„i)=0.30. 
Initial  Yield  Stress  (lo)=  55  MPa,  and  Post  Yield  hardening  law:  0-5, „  =  Yo  +  2.08ePg„ 

Alumina  particle  (Elastic):  Young’s  Modulus  (Ec)=  344.5GPa,  Poisson  Ratio  (1^0)=  0.26,  Critical 
particle  cracking  stress  ((T‘^)=0.3GPa 

The  load  is  applied  in  20  equal  increments  upto  a  total  displacement  of  1  mm  as  shown  in  figure 
4.15a.  For  the  composite  microstructures,  figure  4.16  depicts  the  evolved  macroscopic  and  level- 
2  computational  domains  at  the  end  of  loading  stages.  The  level-0,  level-1,  level-2  elements  are 
shown  in  white,  grey  and  black  respectively  for  the  adapted  macroscopic  mesh  in  4.16a  and  b. 
The  evolution  of  the  levels  and  mesh  with  h-adaptation  is  provided  in  table  3.  The  effect  of  the 
microstructure  size  becomes  more  pronounced  with  increasing  deformation.  A  larger  level-2  domain 
(29  macroscopic  elements)  with  a  smaller  level-0  domain  is  evidenced  for  the  smaller  RVE(i).  The 
effect  of  RVE  size  on  the  pattern  of  particle  cracking  is  very  significant.  The  level-2  region  shows 
24  cracked  particles  for  the  RVE(i)  as  opposed  to  6  cracked  particles  for  the  RVE(ii).  The  path  of 
evolution  of  cracked  particles  is  qmte  different  for  the  two  models.  For  the  RVE(ii),  the  aggregate 
microscopic  cracks  is  found  to  propagate  perpendicular  to  the  macroscopic  loading  direction  and 
all  the  microcracks  have  the  same  orientation  as  the  chain  or  macrocrack.  For  the  smaller  particles 
in  RVE(i),  the  chain  of  microcracks  or  the  macrocrack  is  observed  to  move  at  45°  to  the  loading 
direction  with  individual  microcracks  predominantly  at  90°  to  the  loading  direction. 

The  contour  plots  of  macroscopic  and  microscopic  plastic  strain  distribution  at  the  final  loading 
stage  are  shown  in  figures  4.17  and  4.18  for  the  two  microstructures.  The  macroscopic  strain  distri¬ 
bution  for  both  models  shows  higher  strain  concentration  at  the  corner  of  square  hole  with  increased 
loading.  While  the  macroscopic  strains  are  not  very  different  for  the  two  RVE’s,  significantly  larger 
local  plastic  strains  are  seen  in  the  level-2  microstructure  of  RVE(i).  A  better  representation  of 
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Figure  4.14:  (a)  Macroscopic  asdal  stress-strain  plots  for  a  single  RVE  in  tension  modeled  by  level- 1 
and  level-2  elements;  Maximum  principal  stress  distribution  in  the  inclusion  by  (b)  level-2  element 
and  (c)  level-1  element. 
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Figure  4.16:  Macroscopic  three  level  evolved  me 
loading  cycle  for  (a)  the  larger  RVE  (ii)  and  (b)  ! 
with  the  (c)  the  larger  RVE  (ii)  and  (d)  smaller  I 
with  grey  elements  and  level-2  is  with  black  elem 
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Figure  4.17:  Contour  plot  of  plastic  strain  for  the  composite  square  plate:  (a)  the  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  larger  RVE  (ii) 
model. 


this  difference  is  seen  in  the  graphs  of  figures  4.19  and  4.20.  The  macroscopic  (averaged)  stress  T,s:x 
history  at  the  critical  corner  in  figure  4.19  does  not  indicate  significant  difference  and  hence  exhibits 
little  scale  effect.  The  stress  drops  in  this  figure  correspond  to  particle  damage.  The  histogram 
of  the  evolution  of  particle  cracking  however  shows  a  considerably  different  pattern  and  a  much 
higher  rate  of  cracking  for  RVE(i). 

The  same  problem  is  solved  for  porous  microstructures  with  the  two  sized  RVE’s,  but  with  a 
total  displacement  of  0.4mm.  The  evolved  levels  and  meshes  in  the  computational  models  at  the 
end  of  loading  are  shown  in  figures  4.21a  and  b  and  the  corresponding  level-2  microstructures  in 
figures  4.21c  and  d.  The  evolution  of  the  computational  domain  is  also  charted  in  table  4.  It 
is  interesting  to  note  that  the  difference  in  response  for  the  two  RVE’s  is  insignificant  this  this 
case.  This  may  be  attributed  to  lack  of  matrix  damage  or  softening  in  the  model.  The  plastically 
hardening  matrix  does  not  trigger  adequate  difference  in  the  adaptation  criteria  as  the  particle 
cracking  does  for  composites.  A  larger  level-1  domain  opens  up  with  the  adaptation  criteria  for 
elements  which  appear  to  be  on  the  verge  of  strain  localization,  but  subsequently  do  not  make  the 
transition  to  level-2.  The  contour  plots  of  macroscopic  and  microscopic  plastic  strain  distribution 
in  figures  4.22  and  4.23,  show  similar  macroscopic  strain  distributions,  but  higher  strains  for  the 
microscopic  model  with  smaller  porosity  RVE(i),  due  to  the  proximity  of  voids.  This  again  shows 
the  scale  effect  on  the  solution. 

4.10.2  Homogenization  vs,  multi-level  simulation 

The  effect  introducing  level-2  elements  on  both  macroscopic  and  microscopic  response  is  studied 
by  comparing  a  pure  level- 1  simulation  of  the  square  composite  plate  in  the  previous  example  with 
a  multi-level  simulation.  The  results  shown  are  for  the  larger  particles  in  RVE(ii).  The  figure  4.24 
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Figure  4.20:  Histogram  of  the  number  of  damaged  particles  as  a  function  of  straining. 

shows  the  microstructure  near  the  inside  corner  by  the  two  models  at  the  end  of  loading.  The  boxed 
RVE’s  in  figure  4.24a  symbolize  their  periodic  repetition.  The  periodicity  constraint  results  in  a 
considerably  large  portion  of  the  microstructure  being  damaged  for  the  homogenized  simulation. 
The  direction  of  the  damage  evolution  indicated  by  the  homogenized  model  is  also  different  from 
the  level-2  simulations.  The  stress  along  the  section  A-B  is  plotted  in  figure  4.25  to  evaluate 
the  effect  of  homogenization  on  stress  concentrations  near  the  corner  and  free  edge.  The  two 
models  behave  sinailarly  upto  the  neighborhood  of  the  corner.  While  the  multi-level  model  predicts 
a  higher  peak  near  the  corner,  it  subsides  considerably  to  meet  the  traction  free  edge  conditions. 
The  corresponding  microscopic  level-2  stress  variations  for  the  multi-level  model  are  shown  in  the 
inset. 

4.10.3  Effect  of  heterogeneity  distribution  and  shape 

To  illustrate  the  influence  of  particle  distribution  on  the  macro-microscopic  damage  response,  two 
RVEs  are  selected  with  same  volume  fraction  (20%),  size  (1.0mm)  and  number  of  particles  (25). 
One  has  a  hardcore  distribution  (figure  4.15b  (iii)),  which  is  a  random  distribution  with  a  minimum 
permissible  distance  between  particles,  while  the  other  has  one  cluster  in  a  hardcore  background. 
Proximity  of  particles  within  the  cluster  is  known  to  initiate  damage  in  discrete  microstructures. 
The  starting  macroscopic  mesh  is  the  same  as  in  the  previous  examples  and  a  total  displacement 
of  0.55mm  is  applied  on  the  edge  in  equal  increments.  The  smallest  allowed  size  of  macroscopic 
elements  by  h-adaptation  is  set  to  1mm  such  that  each  E12  element  consists  of  one  RVE. 

The  evolved  macroscopic  models  for  the  two  RVE’s,  at  the  finish  of  the  loading  cycle,  are  shown 
in  figures  4.26a  and  b.  The  level-2  region  (black)  for  the  clustered  RVE  is  larger  than  that  for 
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Figure  4.21:  Macroscopic  three  level  evolved  mesh  for  the  porous  material  at  the  end  of  the  loading 
cycle,  for  (a)  the  larger  RVE  (ii)  and  (b)  smaller  RVE  (i);  The  corresponding  level-2  region  with 
the  (c)  the  larger  RVE  (ii)  and  (d)  smaller  RVE  (i).  Level-0  is  with  white  elements,  level- 1  is  with 
grey  elements  and  level-2  is  with  black  elements. 


Figure  4.25:  Comparison  of  macroscopic  stress  along  section  A-B  by  level-1  and  multi-level 


the  hardcore  RVE.  Within  the  E12  elements,  only  one  element  for  the  hardcore  distribution  expe¬ 
riences  particle  damage  as  shown  in  figure  4.26d.  However,  several  E12  elements  for  the  clustered 
microstructure  exhibit  particle  damage,  mainly  within  the  cluster  (figure  4.26e).  While  the  macro¬ 
scopic  averaged  plastic  strains  show  very  little  difference  for  the  different  microstructures  in  figures 
4.27a  and  4.28a,  the  microscopic  plots  in  figures  4.27b  and  4.28b  clearly  depict  the  influence  of 
distribution.  Much  higher  levels  of  effective  plastic  strain  values  are  observed  within  the  cluster, 
compared  to  significantly  lower  levels  in  the  hardcore  RVE.  Figure  4.30  shows  the  and  evolution 
of  macroscopic  stress  'Sxx  at  the  corner  of  the  square  hole  and  the  number  of  damaged  particles 
as  a  function  of  straining.  The  stress  drops  to  lower  values  for  the  clustered  RVE  due  to  a  larger 
damage  microstructure.  More  than  twice  the  number  of  particles  are  damaged  for  the  clustered 
case  as  shown  in  the  histogram. 

To  investigate  the  influence  of  shape,  a  RVE  (see  figure  4.15b  (iv))  with  the  same  volume  fraction 
(20%),  size  (1.0mm)  and  number  (25)  is  considered.  The  particles  are  elliptical  with  aspect  ratio 
3.0  and  randomly  distributed  and  oriented.  The  evolved  macroscopic  model  in  figures  4.26c  shows  a 
larger  level-2  region  compared  with  other  two  microstructures,  with  several  E12  elements  exhibiting 
particle  damage.  The  much  larger  number  of  cracked  particles  is  also  observed  from  the  histogram 
of  figure  4.30b.  For  this  case,  both  macroscopic  and  microscopic  plastic  strains  in  figure  4.29  are 
also  considerably  larger.  The  macroscopic  stress  shows  a  larger  drop  due  to  the  increased 
damage  in  the  microstructure. 

4.10.4  A  heterogeneous  plate  with  a  macroscopic  holes 

A  different  macroscopic  domain,  viz.  a  plate  with  periodically  repeated  square  diagonal  array  of 
circular  holes  is  considered  in  this  final  example.  The  plate  is  incrementally  loaded  using  prescribed 
displacement  on  the  top  and  bottom  surfaces  to  a  total  extension  of  0.15  mm.  Due  to  periodicity 
and  symmetry,  only  a  part  of  domain  is  modeled  as  shown  in  figure  4.31.  The  radius  of  the  circular 
holes  are  50mm  for  the  100  mm  x  100  mm  square  plate  as  shown  in  figure  4.31.  The  microstructural 
RVE  is  a  20%  volume  fraction  0.4  mm  x  0.4  mm  square  region  with  a  single  circular  particle.  The 
same  material  properties  as  in  the  previous  example  are  used  with  the  only  exception  being  that 
the  critical  paxticle  cracking  stress  cr"'  =  0.2  GPa. 

The  adapted  multilevel  computational  domain  is  shown  in  figure  4.32  and  the  number  of  el¬ 
ements  in  each  levels  with  increasing  are  tabulated  in  table  5.  The  level-2  elements  are  created 
along  a  clear  path  connecting  the  holes  due  to  localization  by  particle  cracking.  Extended  portions 
of  the  damaged  microstructure  in  level-2  regions  are  shown  in  figure  4.32.  The  macroscopic  plastic 
strain  contour  in  figure  4.33a  gives  an  indication  of  ‘hot  spots’  of  evolving  damage  near  the  central 
region.  The  microscopic  strain  plots  in  figure  4.33b  show  a  large  fraction  of  particles  cracked  and 
may  be  interpreted  as  the  initiation  of  localization  to  cause  failure  between  the  holes. 

4.11  Conclusions 

In  this  work,  an  adaptive  multi-level  computational  model  is  proposed  to  (a)  concurrently  predict 
the  evolution  of  stresses  and  strains  at  the  structural  and  microstructural  scales  and  (b)  to  track 
the  incidence  and  propagation  of  microstructural  damage.  The  microstructural  analysis  is  con¬ 
ducted  with  the  Voronoi  cell  finite  element  model  (VCFEM)  for  elastic-plastic  constituents  with 
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Figure  4.26:  Macroscopic  three  level  evolved  mesh  for  the  composite  material  at  the  end  of  the 
loading  cycle  for  (a)  hardcore  distribution  with  circular  particles,  (b)  clustered  distribution  with 
circular  particles,  and  (c)  hardcore  distribution  with  elliptical  particles;  the  corresponding  level-2 
microstructures  for  (d)  hardcore  (e)  clustered  and  (f)  elliptical  RVE’s. 
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Figure  4.29:  Contour  plot  of  plastic  strain  for  the  composite  square  plate:  (a)  the  macroscopic 
averaged  strain  and  (b)  the  level-2  microscopic  strain  at  the  critical  region  for  the  elliptical  RVE 
model. 


Figure  4.30:  The  evolution  of  (a)  at  the  corner  node  of  square  hole  and  (b)  the  number  of 
damaged  particles  for  the  three  microstructures. 
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Figure  4.31:  (a)  Finite  element  model  for  a  part  of  the  composite  plate  with  two  circular  holes  and 
(b)  the  VCFE  model  of  the  microstructural  composite  RVE. 


Figure  4.32:  Macroscopic  three  level  evolved  mesh  for  the  composite  plate  with  two  holes  at  the 
end  of  loading.  The  level-2  regions  in  the  microstructure  are  shown  at  the  two  regions  A  and  B. 
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Figure  4.33:  Contour  plot  of  (a)  the  macroscopic  plastic  strain  and  (b)  level-2  microscopic  plastic 
strain  at  the  region  C  of  the  microstructure 
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particle  cracking.  VCFEM  allows  continuous  changing  element  topology  due  to  progressive  particle 
cracking  with  a  high  level  of  efficiency  and  accuracy  [23,  28].  This  allows  modeling  large  portions 
of  the  complex  microstructure.  A  conventional  displacement  based  elastic-plastic  FEM  code  exe¬ 
cutes  the  macroscopic  analysis.  Adaptive  methods  and  mesh  refinement  strategies  are  developed 
to  create  a  hierarchy  of  computational  sub-domains  with  varying  resolution.  Such  hierarchy  allows 
for  differentiation  between  non-criticaJ  and  critical  regions  and  help  in  increasing  the  efficiency  of 
computations  through  preferential  ‘zoom  in’  regions.  Coupling  between  the  scales  for  regions  with 
periodic  microstructure  is  accomplished  through  asymptotic  homogenization,  while  other  regions 
of  nonuniformity  and  non-periodicity  are  modeled  by  true  microstructural  VCFEM  analysis. 

In  the  hierarchical  model,  three  distinct  levels  evolve  with  progressive  deformation.  In  level- 
0,  a  piecewise  continuous  elastic-plastic  constitutive  law  is  developed  from  a  numerically  created 
database  for  macroscopic  simulations.  The  model  incorporates  exact  microstructural  morphology 
in  its  development  and  hence  is  accurate  from  a  micromechanics  point  of  view.  This  results  in  high 
efficiency  compared  with  two-scale  analysis  by  homogenization.  Level-1  elements,  which  manifest 
with  imminent  damage,  use  asymptotic  homogenization  for  predicting  macroscopic  variables  as 
well  a  variables  in  the  microstructural  RVE.  Simultaneously,  element  refinement  by  h-adaptation 
enables  zooming  in  at  critical  regions  of  intense  deformation.  With  the  rise  in  local  gradients  of 
macroscopic  variables  and  microstructural  damage,  the  assumptions  of  representative  volume  ele¬ 
ment  in  the  microstructure  no  longer  hold.  This  necessitates  a  shift  in  macroscopic  simulations  to 
completely  microscopic  simulations  in  these  regions.  This  region  is  termed  as  the  level-2  where  an 
extended  portion  of  the  microstructure  is  directly  modeled  by  VCFEM.  The  microstructural  model 
is  directly  interfaced  with  level-0  or  level-1  elements  and  a  coupled  analysis  is  performed. 

Several  numerical  examples  are  conducted  with  the  multi-level  model  to  examine  the  effect  of 
various  microstructural  morphologies  on  the  multi-scale  response  of  composite  and  porous  struc¬ 
tural  components.  Specifically  scale  effects,  effects  of  microstructural  distribution  and  shapes, 
and  structural  geometries,  on  the  mechanical  and  damage  response  are  investigated.  The  model 
performs  satisfactorily  in  identifying  critical  regions  and  delineating  the  necessity  of  true  multi¬ 
scale  simulations.  In  conclusions,  the  multi-level  model  developed  reserves  the  potential  to  be  a 
very  effective  tool  in  the  prediction  of  structural  failure  which  can  be  modeled  as  an  aggregate  of 
microscopic  damage. 
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Chapter  5 


Experimental-Computational 
Investigation  Of  Damage  Evolution 
In  Discontinuously  Reinforced 
Aluminum  Matrix  Composite 


5.1  Introduction 

The  commercial  use  of  particle-reinforced  metal  matrix  composites  in  automotive,  aerospace  and 
other  engineering  systems  has  increased  in  the  last  few  decades  due  to  their  potentially  superior 
mechanical  properties,  as  well  as  their  ability  to  reduce  life-cycle  costs  through  enhanced  thermal 
stability  and  weight  reduction.  The  property  advantages  of  these  materials  are,  however,  often  di¬ 
minished  by  general  degradation  of  failure  properties  like  ductility  and  fracture  toughness.  Various 
experimental  and  numerical  studies  [1,  2,  3,  4,  5,  6,  7,  8,  9,  10]  have  been  conducted  to  understand 
the  influence  of  morphological  factors  such  as  volume  fraction,  size,  shape  and  spatial  distribution 
as  well  as  constituent  material  and  interface  properties  on  the  deformation  and  damage  behavior. 
These  studies  have  concluded  that  failure  mechanisms  are  highly  sensitive  to  local  reinforcement 
distribution,  morphology,  size,  interfacial  strength  etc. 

Traditionally  unit  cell  models  [11,  12,  13,  14,  15]  based  on  the  finite  element  analysis  have 
been  used  to  predict  the  onset  and  growth  of  evolving  damage  in  composite  materials.  While  these 
models  provide  valuable  insights  into  the  microstructural  damage  processes,  simple  morphologies 
idealize  actual  microstructures  for  many  engineering  materials  that  bear  little  relationship  to  the 
actual  stereographic  features.  These  deficiencies  have  been  circumvented  in  [16,  8,  17],  where  com¬ 
putational  models  of  discontinuously  reinforced  materials  with  random  spatial  dispersion  have  been 
considered.  Richmond  and  coworkers  [18, 19]  have  investigated  the  effect  of  morphology  on  damage 
in  composite,  porous  and  polycrystalline  materials  by  modeling  actual  geometries  obtained  from 
2D  micrographs.  Using  the  Voronoi  Cell  finite  element  model,  Ghosh  et.  al.  [9,  10]  have  examined 
the  effect  of  various  spatial  dispersions  and  particle  shape  and  size  on  the  damage  initiation  and 
evolution  process  in  ductile  matrix  composites. 

Many  characterization  studies  with  2-D  microstructures  e.g.  [20,  21,  22,  23]  have  also  been 
conducted  to  understand  the  relation  between  microstructural  morphology  and  damage.  Experi- 
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mental  research  [25,  26,  27,  28,  29]  has  however  pointed  to  the  necessity  of  examining  the  full  3D 
characteristics  for  understanding  the  damage  process.  These  studies  infer  that  2-D  assessment  can 
sometimes  be  misleading,  especially  in  the  presence  of  spatial  clustering.  Non-destructive  evalua¬ 
tion  methods,  e.g.  techniques  based  on  ultrasonics  e.g.  [30],  acoustic  emission  [4]  and  X-ray  based 
computer  tomography  (CT)  [31,  32]  have  emerged  as  potential  methods  for  studying  3D  damage. 
However  many  of  these  systems  are  thus  far  not  capable  of  achieving  spatial  resolutions  required  to 
accurately  capture  microscopic  particles  and  damage  in  the  particle  reinforced  MMC’s.  Buffiere  et. 
al  [33]  are  developing  a  CT  technology  to  yield  tomographic  images  with  a  higher  spatial  resolution. 

This  work  deals  with  a  combined  experimental-computational  approach  to  study  the  evolution 
of  microscopic  damage  to  cause  complete  material  failure  in  commercial  SiC  particle  reinforced 
aluminum  alloys  or  DRA’s.  Through  a  combination  of  2D  and  3D  characterization  and  analysis 
models,  it  is  intended  to  understand  what  aspects  of  microstructural  morphology  that  are  most 
critical  for  damage  nucleation  and  evolution.  Since  it  is  difficult  to  identify  the  microcrack  growth 
process  once  a  material  has  failed  completely,  an  interrupted  testing  technique  is  designed.  Sub¬ 
sequently,  sample  microstructures  in  the  severely  necked  region  are  microscopically  examined  in 
3D  using  a  serial  sectioning  method  discussed  in  [24,  25,  26].  Computer  simulated  equivalent  mi¬ 
crostructures  are  tessellated  into  meshes  of  2-D  and  3-D  Voronoi  cells.  Various  characterization 
functions  of  geometric  parameters  are  generated  and  a  sensitivity  analysis  is  conducted  to  explore 
the  influence  of  morphological  parameters  on  damage.  2D  characterization  functions  are  compared 
with  3D  to  evaluate  the  effectiveness  of  modeling  the  2D  micrographs.  Modeling  of  the  initia¬ 
tion  and  propagation  of  damage  is  conducted  with  Voronoi  Cell  Finite  Element  Method  (VCFEM) 
[9,  10,  34,  35].  Each  Voronoi  cell  element  may  consist  of  a  matrix  phase,  an  inclusion  phase  and 
a  crack  phase.  Damage  initiation  by  particle  cracking  is  assumed  to  follow  a  maximum  principal 
stress  based  Rankine  criterion.  The  VCFEM  for  particle  cracking  has  shown  a  significant  promise 
in  modeling  large  aggregates  of  heterogeneities.  While  the  appropriateness  of  3D  analyses  is  rec¬ 
ognized  for  this  study,  the  3D  VCFEM  (under  development)  does  not  currently  have  all  necessary 
features.  Due  to  enormous  computing  requirements  of  conventional  3D  FEM  models,  various  stud¬ 
ies  have  resorted  to  simplified  manifestations  of  complex  geometries  and  properties  e.g.  [7,  41,  15]. 
This  study  is  restricted  to  2D  in  the  form  of  VCFEM  analyses  of  section  micrographs.  Finally,  the 
effect  of  size  and  characteristic  lengths  of  representative  material  element  (RME)  on  the  extent  of 
damage  in  the  model  systems  is  also  investigated. 


5.2  Experiments  for  Damage  Assessment 

5,2.1  Interrupted  tests 

The  material  analyzed  in  this  work  is  a  discretely  reinforced  commercial  aluminum  that  is  fabricated 
by  a  powder  metallurgy  process  [36].  It  consists  of  extruded  commercial  X2080  aluminum  alloy 
with  15%  volume  fraction  SiC  particles.  The  X2080  matrix  has  a  nominal  alloy  composition  with 
weight  percentages  of  3.8%  Cu,  1.8%  Mg  and  0.2%  Zr,  in  addition  to  low  impurity  contents  of  Fe 
and  Si.  The  precipitation  hardened  X2080  aluminum  alloy  system  is  naturally  aged  by  heat  treating 
for  4  hours  at  930^^,  followed  by  cold  water  quench  and  aging  for  2  days  at  room  temperature. 

An  important  object  in  this  failure  study  is  to  obtain  adequate  microstructural  data  that  depict 
the  growth  of  damage  into  a  major  failure  path.  In  general,  it  is  difficult  to  identify  the  domi¬ 
nant  damage  mechanisms  and  also  the  microcrack  growth  process,  once  a  material  has  fractured 
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completely.  Thus  an  interrupted  testing  technique  is  designed  where  the  load  and  deformation  are 
halted  in  the  material  instability  zone,  following  necking  but  prior  to  fracture.  The  tests  assume 
that  the  major  cracks  are  essentially  prominent  is  this  stage,  and  are  helpful  in  understanding  the 
linkage  mechanism  of  microcracks  or  particle  debonds  to  facilitate  growth  of  the  dominant  damage. 
To  initialize  the  testing,  estimates  of  the  necking  and  fracture  strains  are  first  obtained  by  observing 
the  behavior  of  a  tension  test  to  failure.  The  uniaxial  tension  tests  are  executed  on  a  MTS  810  ma¬ 
terial  system  with  a  HP  7044  X-Y  recorder  to  monitor  the  loads  and  strains,  and  the  critical  strains 
are  measured  with  a  MTS  632.11  strain  gauge  extensometer.  Following  the  initialization,  strain 
controlled  interrupted  tests  are  carried  out,  in  which  the  specimens  are  loaded  to  the  instability 
region  before  the  load  is  stopped. 

Figure  5.1a  shows  a  typical  tension  specimen  for  the  naturally  aged  material.  Data  for  six 
specimens  of  this  material,  viz.  tl,  t2,  t3,  t4,  t5  and  t6  are  tabulated  in  table  1.  The  specimens 
tl,  t2,  t4  and  t6  are  obtained  from  the  outer  annulus  region  of  the  stock  material  while  t3  and 
t5  are  from  the  central  core  regions.  The  initialization  of  the  test  to  study  the  entire  material 
behavior  and  estimate  the  post-instability  region  is  done  with  specimens  tl  and  t2.  The  material 
load- displacement  curve  is  plotted  in  figure  5.1b,  from  which  the  necking  strain  is  obtained  from 
the  peak  load  value.  For  the  specimen  tl,  the  test  is  conducted  at  a  strain  rate  is  e  =  5xl0“^sec“^ 
and  the  necking  strain  and  fracture  strain  are  found  to  be  =  9.15%  and  e,-  =  9.40%  respectively. 
The  short  instability  region  in  tl  prompts  a  reduced  strain  rate  e  =  3xl0“‘*sec“^  for  specimen  t2, 
for  which  =  9.05%  and  e,-  =  9.20%. 

In  table  1,  e,  €„  and  e,-  correspond  to  the  strain  rate,  the  necking  strain  and  the  interrupted 
strain  respectively.  The  interrupted  strain  coincides  with  the  fracture  strain  in  the  event  that 
fracture  precedes  the  load  stoppage.  This  is  indicated  with  F  or  I  in  the  table.  Load  interruption  is 
only  possible  for  the  specimens  t3  and  t6  due  to  the  extremely  short  post-instability  range  of  this 
material  in  comparison  with  the  resolution  of  the  loading  mechanism.  The  necking  strains  for  the 
specimens  tl,  t2,  t4  and  t6  are  in  the  range  of  9.00%  ~  9.30%,  while  those  for  specimens  t3  and  t5 
are  in  the  9.80%  ~  10.20%  range.  This  difference  is  possibly  due  to  gradients  created  by  the  heat 
treatment  at  different  locations  in  the  stock  material.  The  core  cools  slower  and  more  uniformly 
regions  near  the  surface.  This  results  in  the  more  uniform  microstructure  and  larger  necking  and 
fracture  strains  for  specimens  (t3  and  t5)  located  near  the  core  of  the  stock  material. 

5.2.2  Damage  examination  and  microscopic  analysis 

To  examine  the  dependence  of  microstructural  damage  on  the  local  morphology,  serial  sectioning 
of  sample  coupons  extracted  from  the  load-interrupted  specimens  t3  and  t6,  is  invoked.  This 
method,  discussed  in  [24,  25,  26],  involves  gradual  removal  of  material  layers  to  obtain  a  series 
of  scanning  electron/optical  micrographs,  representing  sections  of  a  microstructure.  It  is  a  very 
effective  method  for  reconstructing  3D  microstructures  from  a  series  of  2D  sections  of  particulate 
reinforced  composites,  requiring  a  resolution  of  few  microns.  Prior  to  sectioning,  locations  are 
selected  in  figure  5.1a  for  cutting  out  the  sample  coupons.  X-rays  and  acoustic  microscopy  with 
a  AEROTECH  UNIDEX  11  acoustic  microscope  with  a  resolution  of  about  50/im  are  used  in 
this  process  to  detect  regions  that  contain  major  crack  paths.  Polished  surfaces  of  these  extracted 
samples  are  then  examined  by  a  Nikon  optical  microscope  for  major  damage  sites.  For  the  specimen 
t3,  shorter  cracks  passing  through  2  ~  3  particles  at  most  are  found.  However,  for  the  specimen 
t6,  a  larger  crack  passing  through  5  ~  6  particles  is  identified,  and  is  consequently  chosen  for 
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Figure  5.1;  (a)  Interrupted  uniaxial  tensile  test  specimen  for  naturally  aged  material  and  sample 
coupon  for  serial  sectioning(unit  in  mm);  (b)  Load-strain  plots  for  two  specimens  of  naturally  aged 
DRA.  Dark  points  indicate  where  the  loading  is  interrupted  or  where  the  specimen  is  fractured. 
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analysis.  Coupons  of  approximate  size  6  mm  x  6  mm  x  6  mm  are  subsequently  prepared  for 
the  serial  sectioning  operation  to  sequentially  expose  parallel  sections  of  the  microstructure.  As 
discussed  in  [24,  25,  26],  parallel  layers  in  a  direction  perpendicular  to  the  straining  direction  (see 
figure  5.1a)  are  removed  using  a  precision  dimple  grinder.  The  depth  of  material  removal  per  step 
is  selected  such  that  each  particle  is  sectioned  at  least  once,  ensuring  that  all  particles  of  interest 
are  adequately  captured  in  the  micrographs.  For  the  DRA  considered,  the  particle  size  range  is 
approximately  3-25  /im,  with  an  average  size  of  ~  9.2  fim  and  the  standard  deviation  is  3.891  fim. 
The  section  to  section  step  size  is  chosen  to  be  2  /xm,  corresponding  to  a  total  traversed  thickness 
of  36/im  for  18  sections.  Two  typical  micrographs  showing  damage  are  depicted  in  figures  5.2, 
for  which  the  horizontal  corresponds  to  the  loading  direction.  The  micrographs  are  then  serially 
stacked  using  a  graphic  software  [37]  to  yield  3D  microstructures  as  shown  in  figure  5.3a.  The 
precise  3D  location,  shape,  size  and  orientation  of  each  particle  can  be  obtained  at  a  fairly  high 
resolution  by  this  method. 

5.2.3  Major  observations 

The  micrographs  of  serial  sections  3  and  5  in  figure  5.2,  perpendicular  to  the  middle  plane  of  the 
tensile  specimen,  provide  important  information  on  the  evolution  of  the  dominant  damage  path 
in  the  material.  A  dominant  damage  path  is  clearly  seen  in  the  boxed  regions.  The  damage  size 
progressively  diminishes  with  increasing  sections,  indicating  the  end  of  the  cracked  particles.  The 
particle  area  fraction  (AF),  total  number  of  particles  (NP)  and  total  number  of  cracked  particles 
(NCP)  for  each  section  micrograph  are  presented  in  table  2.  Generally  speaking,  sections  with  large 
AF  and  NCP  are  found  to  contain  the  larger  cracks.  The  3D  image  by  assembling  2D  micrographs 
in  figure  5.3a  also  shows  the  dominant  damage  path  in  the  boxed  region. 

From  the  microscopic  observation  results,  it  is  found  that  for  the  naturally  aged  material,  the 
main  mode  of  damage  is  by  particle  cracking.  Large  particles  in  particle  rich  regions  are  more 
susceptible  to  cracking  than  those  in  particle  sparse  regions.  Microcracks  in  the  particle  rich  areas 
link  up  to  form  paths  of  dominant  damage.  The  linkage  and  evolution  of  these  larger  cracks  lead 
to  the  overall  failure  of  the  material.  These  paths  are  approximately  perpendicular  to  the  tensile 
loading  direction.  Thus,  spatial  distribution  of  particles  plays  a  more  important  role  in  damage 
than  particle  size  for  this  material. 


5.3  Equivalent  Microstructure  &  Mesh  Generation 

The  actual  3D  geometry  of  particles,  as  seen  in  figure  5.3a,  can  be  quite  complex  and  an  exhaustive 
database  is  required  to  store  all  geometric  details.  To  avert  this,  equivalent  microstructures  that 
closely  approximate  the  actual  morphology  but  are  computationally  less  demanding,  are  generated. 
In  this  process,  each  particle  and  an  associated  crack  are  replaced  by  equivalent  ellipses  (in  2D)  or 
ellipsoids  (in  3D).  This  method  economizes  the  image  analysis  and  characterization  process  by  way 
of  well  known  geometric  properties.  For  obtaining  equivalent  microstructures,  digitized  image  data 
is  first  transferred  into  a  binary  format  to  distinguish  between  the  particle,  matrix  and  crack  phases. 
The  0th  (Jo),  1st  and  2nd  (Ixx^Iyy)  order  geometric  moments  are  then  computed  for  each 

particle  by  adding  contributions  from  each  voxel  (in  3D)  or  pixel  (in  2D)  that  lies  within  the  particle 
boundary.  For  2D  microstructures  the  computed  moments  are  equated  to  the  moment  formulae 
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Figure  5.2:  Micrographs  of  diflFerent  sections  of  the  t6  specimen  showing  cracked  particles;  (a) 
section  3,  (b)  section  5. 
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for  ellipses  to  evaluate  the  centroidal  coordinates  (a:c,2/c)»  half  major  and  minor  axis  lengths  (a,  6), 
and  orientation  6  of  the  major  axis  from: 


Xc 

e 


9.  =  |.  »  =  \/5(Ci  +  \/cr^,  J=\/5(c,-\/cJ-4C2) 
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(5.1) 


where  Ci  =  4*  ( — ^75^)  and  C2  =  For  3D  microstructures,  the  centroidal  coordinates 
(^cii/ct^c)  of  the  equivalent  ellipsoid  are  first  evaluated  from  the  0th  and  order  moments  as: 

^  j  yc  =  ^  ~  ^  principal  directions  (or  orientations  of  the  three  axes)  for  the 

ellipsoids  are  obtained  from  the  eigen-values  of  the  2nd  order  moments  Ijj  {i  =  1..3,  j  =  1..3).  The 
major  (2a),  intermediate  (26)  and  minor  (2c)  axes  of  the  equivalent  ellipsoids  are  then  obtained 
from  the  principal  moments  hyhfh  as: 


=  +  C=  J|(/i  +  /2-/3)  (5.2) 

A  simulated  3D  microstructure  with  particles  (grey)  and  cracks  (black)  is  shown  in  figure  5.3b. 
The  microstructures  are  then  tessellated  into  a  mesh  of  2D  and  3D  Voronoi  cells,  by  surface  based 
algorithms  detailed  in  [25,  26].  In  figure  5.3c,  the  mesh  of  Voronoi  cells  is  created  based  on  the 
morphology  of  particles,  while  in  figure  5.3d  the  mesh  is  due  to  tessellation  based  on  the  geometry 
of  particle  cracks.  Tessellation  into  a  mesh  of  Voronoi  cells  plays  an  important  role  in  developing 
geometric  descriptors  for  quantitative  characterization.  They  represent  regions  of  immediate  influ¬ 
ence  for  each  heterogeneity  and  also  define  neighbor  of  each  heterogeneity  from  individual  faces  or 
edges  of  the  Voronoi  cells.  This  facilitates  easy  evaluation  of  parameters  like  local  area  fractions, 
near  neighbor  and  nearest  neighbor  distances  and  orientations. 


5.4  Microstructure  and  Damage  Characterization 

The  morphology  of  particles  and  associated  damage  or  microcracks  can  be  characterized  by  various 
functions  of  size,  shape,  orientation  and  spatial  distribution.  A  number  of  these  classifier  functions 
have  been  used  by  the  authors  and  others  in  [24,  25,  26,  34,  35,  20,  21,  22,  6]  to  characterize  various 
aspects  of  microstructural  morphology.  In  this  section,  some  of  these  functions  are  considered 
for  the  3D  microstructure  and  2D  micrographs  to  investigate  the  relation  between  morphological 
characteristics  and  the  path  of  dominant  damage  in  the  material.  The  specimen  t6  with  a  large 
microcrack  is  considered  for  this  study. 

In  the  first  exercise,  a  sensitivity  analysis  is  done  with  the  simulated  3D  microstructure  in  figure 
5.3b  to  reveal  the  dependence  of  damage  on  microstructural  variables.  Two  damage  parameters, 
viz.  the  number  fraction  of  cracked  particles  (nf)  and  the  volume  fraction  of  cracked  particles  (vf) 
are  chosen  to  manifest  the  damage  level  in  the  DRA.  Six  microstructural  parameters  are  considered, 
viz.  (i)  particle  equivalent  size  (diameter);  (ii)  nearest  neighbor  distance  computed  as  the  distance 
between  particles  that  share  a  common  Voronoi  cell  edge;  (iii)  local  volume  fraction  measured  as 
a  ratio  of  the  particle  size  to  that  of  the  associated  Voronoi  cell;  (iv)  particle  shape  or  ellipsoid 
aspect  ratio;  (v)  nearest  neighbor  orientation,  measured  as  the  angle  between  a  line  joining  the 
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centers  of  particle  and  its  nearest  neighbor,  and  the  loading  direction;  and  (vi)  particle  orientation 
with  respect  to  the  loading  direction.  The  cracked  particle  fractions  are  plotted  as  functions  of 
these  parameters  in  figure  5.5.  A  linear  interpolation,  obtained  by  a  least  square  fit,  yields  the 
corresponding  overall  gradient  or  slope. 

While  both  the  nf  and  vf  plots  coincide  for  the  particle  size  plot  (i),  large  differences  are  noted 
for  nearest  neighbor  distances  (ii)  and  aspect  ratios  (iv).  Largest  slopes  of  these  plots  are  observed 
with  particle  size,  nearest  neighbor  distance  and  local  volume  fraction.  This  infers  that  the  strongest 
influence  on  particle  cracking  comes  from  the  size  and  local  spatial  distribution.  Particle  shape  has 
a  relatively  smaller  effect  on  damage  initiation.  Sensitivity  of  damage  to  particle  orientation  and 
nearest  neighbor  orientation  is  found  to  be  minimal  for  this  material. 

The  characteristics  of  particles  forming  the  dominant  damage  path  (within  the  marked  box  in 
5.3a)  are  compared  with  those  for  all  cracked  particles  in  the  histograms  of  figure  5.6.  The  dotted 
lines  correspond  to  all  cracked  particles  while  the  shaded  areas  are  for  cracked  particles  in  the 
dominant  damage  region  only.  The  histograms  are  with  respect  to  three  variables  that  are  found 
to  play  important  roles  in  the  damage  process,  viz.  the  particle  size,  nearest  neighbor  distance 
and  orientation  with  respect  to  the  loading  direction.  While  the  range  of  sizes  for  all  cracked 
particles  is  4  ^  13/im,  that  for  the  particles  forming  the  dominant  damage  path  is  5.7  ^  ISfim, 
This  reveals  that  larger  particles  generally  contribute  to  dominant  damage  path.  The  plot  for 
nearest  neighbor  clearly  exhibits  the  influence  of  particle  rich  areas  (clustering  or  alignment)  on 
the  preferential  growth  of  damage.  The  nearest  neighbor  distance  for  particles  in  the  dominant 
damage  path  are  in  the  range  0.4  ^  3.7/im  when  compared  with  the  range  0.4  12.3//m  for  all 

cracked  particles.  The  histogram  of  cracked  particles  as  a  function  of  the  orientation  with  respect  to 
the  loading  direction  reveals  that  particles  with  major  axis  along  the  loading  direction  (0^  and  180°) 
are  generally  susceptible  to  cracking.  This  is  much  more  pervasive  for  particles  in  the  dominant 
damage  path,  due  to  the  smaller  cross-sectional  areas  normal  to  loading.  In  conclusion,  particles 
in  the  dominant  damage  path  generally  have  larger  size,  are  in  particle  rich  areas,  and  are  oriented 
in  loading  direction. 

Finally,  it  is  of  interest  to  identify  discriminating  characteristics  of  2D  micrographs  that  may 
be  helpful  in  making  dominant  damage  predictions  for  the  actual  3D  microstructures.  Two  repre¬ 
sentative  micrographs,  viz.  section  1  which  contains  a  dominant  damage  and  section  14  with¬ 
out  any  dominant  damage,  but  only  scattered  particle  cracks  are  compared  with  the  3D  mi¬ 
croregion.  Four  important  characterization  functions  viz.  (a)  the  probability  density  function 
of  particle  equivalent  size  (diameter),  (b)  the  probability  density  function  of  the  nearest  neighbor 
distance,  (c)  the  probability  density  function  of  the  local  area/volume  fraction  and  (d)  a  trans¬ 
formation  function  L{r)  of  a  second  order  intensity  function  K{r),  are  plotted  in  figure  5.7  for 
2D  and  3D  micrographs.  The  second  order  intensity  function  K{t)  and  its  transformed  functions 
in  2D,  and  L{r)  =  in  3D)  capture  second  order  statistics  of  spatial 

distributions  are  used  as  a  graphical  tool  for  detecting  departures  from  a  homogeneous  Poisson 
process  [34,  35,  25,  26].  The  plot  of  L{r)  vs.  r  is  a  45°  straight  line  for  a  pure  Poisson  distribution. 

The  plots  distinctly  reveal  a  few  important  features  of  the  micrographs.  The  particle  size 
distribution  for  the  two  2D  micrographs  are  similar  and  the  tails  are  significantly  shorter  than 
3D.  As  is  expected,  3D  particle  sizes  are  larger  than  2D  particle  section  sizes  due  to  sectioning 
along  non-principal  planes.  However  the  probabilities  of  both  the  nearest  neighbor  distances  and 
local  area  fractions  in  figures  5.7b  and  c  yield  a  distinguishing  characteristic.  The  micrograph  with 
dominant  damage  has  peaks  and  valleys,  as  well  as  tails  that  are  very  similar  to  that  for  3D.  The 
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peaks  which  reflect  particle  rich  regions  and  the  tails  which  reflect  sparse  areas  are  both  found  to 
be  important  discriminants.  Deviation  from  the  L{r)  =  r  function  or  the  45°  line  represents  a 
bias  towards  clustering.  The  section  with  the  dominant  damage  has  a  larger  deviation  from  the 
random  distribution  in  comparison  with  the  section  without  major  cracking,  and  is  closer  to  the 
3D  response.  In  summary,  it  may  be  concluded  that  when  analyzing  2D  sections,  the  liklehood  of 
better  representation  of  dominant  damage  are  for  those  sections  that  have  higher  peaks  at  lower 
near  neighbor  distances  with  longer  tails  and  have  higher  deviation  from  the  Poisson  distribution. 
Similar  observations  have  also  been  made  in  [5,  1,  38,  29]. 


5.5  Damage  Simulation  by  Voronoi  Cell  FEM 

Two  dimensional  plane  strain/stress  simulations  of  the  microstructural  damage  evolution  is  con¬ 
ducted  by  the  Voronoi  cell  finite  element  model  (VCFEM)  described  in  [9, 10,  34,  35].  The  current 
2D  VCFEM  only  accommodates  particle  cracking,  and  hence  matrix  cracking  is  ignored  in  the  sim¬ 
ulations.  The  simulations  are  useful  in  understanding  the  damage  evolution  process  by  a  sequence 
of  particle  cracking.  Rectangular  195/rm  x  155.018/im  micrographs  as  shown  in  figure  5.9a,b  are 
analyzed  with  monotonically  increasing  strains.  Periodicity  boundary  conditions  are  imposed  by 
requiring  edges  to  remain  straight  and  parallel  to  the  original  direction  throughout  deformation  as: 

Ux  -  0  {on  X  =  0)  ,  Uy  =  0  (on  y  =  0)  ,  Ux  =  Uap  {on  x  =  Lx)  ,  Uy  =  D*  {on  y  =  Ly) 

Ty  =  Q  {on  X  =  0/Lx)  ,  Tx  =  0  {on  y  =  0/Ly)  (5.3) 

where  Uap  is  an  applied  displacement  and  D*  is  determined  from  the  average  force  condition 
Txdx  =  0  on  y  =  Ly.  The  reinforcing  phase  of  SiC  particles  are  assumed  to  be  brittle  and 
is  modeled  with  the  linear  elastic  properties:  Young’s  modulus  E  =  427  GPa,  Poisson’s  ratio  v 
=  0.17.  The  aluminum  matrix  material  is  assumed  to  be  ductile  and  is  modeled  by  small  defor¬ 
mation  isotropic  hardening  elasto-plasticity  theory  with  properties:  Young’s  modulus  E  =  72 
GPa,  Poisson’s  ratio  i>  =  0.33  and  the  post  yield  elastic-plastic  behavior  is  obtained  from  [39]  as 
shown  in  figure  5.8.  Microstructural  damage  by  particle  cracking  is  assumed  to  be  governed  by 
a  maximum  principal  stress  or  Rankine  criterion.  In  this  criterion,  a  crack  is  initiated  when  the 
maximum  principal  stress  in  tension  exceeds  a  critical  fracture  stress  Ccr  at  a  point.  The  crack  is 
oriented  at  right  angle  to  the  principal  stress  direction.  The  critical  stress  Ocr  is  also  influenced  by 
the  particle  size  due  to  the  existence  of  microcracks.  To  account  for  the  size  effect  in  (Tct,  a  WeibuU 
distribution  based  criterion  is  used,  where  the  probability  of  particle  fracture  Pf{A,  cr)  is  related  to 
the  particle  volume/area  v  and  the  maximum  principal  stress  oj  as: 

Pf{'>’,<^l)  =  1  -  exp[-v{—)'^]  (5.4) 

where  oq  and  m  are  two  material  parameters  in  the  WeibuU  distribution  that  are  calibrated  from 
experiments. 
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5.5.1  Calibration  of  Weibull  parameters  (Tq  and  m 

In  the  two  parameter  Weibull  model,  the  fraction  of  fractured  particles  may  be  obtained  (see 
[41,  42,  26])  from  a  known  probability  distribution  of  particle  volumes  p{v)^  as: 

tVmax  ^  7)* 

p(v)  =  /  p(v)  Pf(v,(Ti)dv  «  (^-5) 

•'Vmin  "0 

where  p(v,)  is  the  probability  density  distribution  of  particle  volume/area  u,-.  The  entire  area  is 
divided  into  N  intervals  such  that  Au,-  =  u,-  —  u,_i,  <t}  is  the  average  particle  maximum  principal 
stress  for  particles  with  size  in  the  range  of  [u,_i ,  u,]  and  vq  is  a  reference  area  taken  to  be  the  average 
area.  The  fraction  of  cracked  particles  p  is  readily  obtained  from  the  experimental  micrographs. 
Again,  the  section  micrographs  2,  8  and  14  are  used  to  calibrate  the  Weibull  parameters.  The 
fractions  of  cracked  particles  and  the  average  particle  area  for  these  three  sections  are  31.78%, 
24.76%,  28.57%  and  53.43,  48.91  and  52.67  pm?  respectively.  The  maximum  principal  stress  a\  for 
each  particle  is  obtained  from  VCFEM  simulation  prior  to  the  onset  of  particle  cracking  at  a  true 
strain  of  c  =  8.88%.  From  the  experimental  observations  it  is  assumed  that  no  major  damage  has 
initiated  at  this  strain.  The  Weibull  parameter  m  is  assumed  to  take  integer  values  between  1  and 
8  following  [42,  26]  and  the  corresponding  values  of  (Jq  are  given  in  table  3. 

The  Weibull  parameters  are  also  calibrated  using  a  3D  ABAQUS  model  simulation  of  a  cubic 
unit  ceU  with  a  single,  15%  volume  fraction,  spherical  particle  as  described  in  [26].  The  1x1x1  unit 
cell  model  has  a  particle  of  radius  R  =  0.66.  A  modified  form  of  equation  5.5  is  used  to  account 
for  the  shape  variability  of  the  particles  as 

t^max  t^max 

p{a,v)=  /  p{a)p{v)  Pf{v,(Ti)dv  (5.6) 

where  a  corresponds  to  the  particle  aspect  ratio.  The  particle  size  and  shape  distribution  functions 
p{v)  and  p(a)  are  calculated  from  the  computer  simulated  representation  of  the  actual  3D  mi¬ 
crostructure  shown  in  figure  5.3.  This  average  particle  volume  and  the  fraction  of  cracked  particle 
are  directly  computed  from  figure  5.3  as  v  =642.0pm^  and  p  =  45.48%.  The  average  particle  stress 
at  a  macroscopic  strain  e  =  8.88%  is  obtained  from  the  ABAQUS  simulation  as  Up  —  862.60  MPa. 
Results  of  calibration  with  and  without  shape  effects  are  documented  in  table  3.  It  is  found  that  the 
best  agreement  in  ctq  for  all  2D  sections  and  3D  is  obtained  for  m  between  4  and  5.  Consequently 
the  parameter  is  chosen  to  be  m  =  4.2.  The  corresponding  value  of  oq  for  section  2  is  3.04  GPo,  for 
section  8  is  3.19  GPa,  for  section  14  is  2.79  GPa  and  the  average  of  these  sections  is  ao  =  S.OlGPa. 

Results  of  VCFEM  analysis  of  the  simulated  micrographs  of  sections  1,  3,  5  and  9  are  provided 
in  table  4.  The  number  of  cracked  particles  at  a  macroscopic  strain  of  8.88%  by  VCFEM  are 
compared  with  experimental  results.  While  the  general  agreement  is  quite  good,  it  is  seen  that  the 
concurrence  is  particularly  favorable  when  the  simulation  is  conducted  with  a  cq  that  is  obtained 
from  a  section  that  is  near  to  the  one  being  analyzed.  For  example,  the  results  of  sections  1  and 
3  are  very  good  when  (Tq  =  3.01  GPa,  which  is  obtained  from  section  2.  This  concurrence  may 
be  attributed  to  the  similarity  in  the  distribution  of  heterogeneities  in  neighboring  sections,  and 
suggests  that  spatial  distribution  has  a  strong  effect  on  the  WeibuU  parameters. 
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Figure  5.5:  Histograms  comparing  characteristics  of  particles  with  dominant  damage  with  all 
cracked  particles 


Microscopic  Damage  Analysis 

Various  results  for  section  1  which  contain  a  dominant  damage  path  are  generated  by  VCFEM 
simulation  and  compared  with  experimental  observations  in  this  section.  The  macroscopic  stress- 
strain  plot  for  plane  strain  and  plane  stress  assumptions  are  compared  with  experimental  results  in 
figure  5.8a.  The  overall  yield  strength  is  better  predicted  by  the  plane  stress  model.  However,  the 
post  yield  behavior  with  plane  strain  conditions  is  much  closer  to  the  experimental  results.  The 
initial  higher  yield  strength  is  expected  with  plane  strain  due  to  the  plastic  constraint  caused  by 
the  =  0  condition.  A  shifted  stress-strain  plot  (modified  plane  strain  VCFEM  result  in  figure 
5.8a)  where  the  stresses  are  reduced  by  the  initial  difference  in  yield  stress  shows  a  very  good 
match  between  experiments  and  simulation.  Thus  plane  strain  assumptions  are  used  in  subsequent 
computations.  Figure  5.8b  is  intended  to  predict  the  onset  of  plastic  instability  by  the  model  and 
compare  it  with  actual  fracture  observed  in  the  experiments.  The  use  of  the  Considere  criterion 
to  predict  the  onset  of  plastic  instability  has  been  suggested  by  Llorca  [2,  41]  in  the  absence  of 
dilatational  strain  associated  with  reinforcement  fracture.  In  this  criterion,  the  average  stress  a  is 
related  to  the  strain  hardening  rate  ^  as 


(5.7) 


The  strain  derived  from  this  relation  corresponds  to  the  lower  bound  of  the  tensile  ductility  since 
it  controls  the  composite  load  bearing  capacity.  Three  sets  of  curves  are  plotted  in  the  figure  5.8b 
corresponding  to  the  matrix  material,  the  VCFEM  results  in  plane  strain  and  the  experimental 
results.  It  is  seen  that  the  Considere  criterion  (junction  of  the  two  curves)  predicts  the  experimental 
point  corresponding  to  the  onset  of  fracture  rather  well.  Additionally  the  2D  prediction  of  the  plane 
strain  simulation  is  also  quite  good  and  can  be  used  with  reasonable  confidence. 

Microstructural  results  of  the  simulation  are  compared  with  experiments  in  figure  5.9.  The 
computed  micrograph  with  evolved  damage  for  section  1  is  compared  with  the  experimental  micro¬ 
graph  at  8.88%  strain  in  figure  5.9a  and  b.  The  damaged  particles  are  shown  with  the  contained 
crack.  Most  damaged  particles  in  the  simulation  coincide  with  the  experimental  results,  with  the 
box  indicating  the  dominant  damage  path.  The  damage  path  is  approximately  perpendicular  to  the 
tensile  loading  direction.  Figure  5.9b  also  shows  the  contour  plot  of  effective  plastic  strain  in  the 
ductile  matrix  and  indicates  the  path  of  damage  linkage.  The  plastic  strain  is  higher  and  localized 
between  cracked  particles  and  this  is  expected  to  cause  matrix  cracking.  The  number  fraction  of 
cracked  particles  in  different  size  ranges  are  plotted  in  figure  5.9c.  Very  good  agreement  is  seen 
between  simulation  and  experimental  results.  Figure  5.9d  is  a  contour  plot  of  the  particle  fracture 
probability  at  8.88%  strain.  The  black  shade  corresponds  to  the  highest  probability  and  fractured 
particles  are  illustrated  in  white  with  a  crack.  Similar  plots  (not  shown)  at  earlier  stages  of  defor¬ 
mation  show  that  several  particles  with  higher  probability  at  the  smaller  strain  have  cracked  with 
deformation.  The  number  fraction  of  cracked  particles  as  a  function  of  straining  are  plotted  for 
sections  1,  5  and  9  together  with  the  experimental  observation  in  figure  5.10.  At  lower  strains  the 
number  fraction  of  of  cracked  particles  for  sections  1  and  5  with  particle  rich  regions  are  higher  than 
that  for  section  9.  This  is  due  to  higher  stress  concentrations  particle  rich  areas  that  are  enough  to 
fracture  some  particles  even  at  low  strains.  With  increasing  strains  more  particles  start  to  crack  in 
the  section  9  and  exceeds  that  for  section  5  which  has  less  particles  in  the  clustered  regions.  The 
2D  simulations  however  exhibit  less  cracked  particles  than  that  in  the  actual  3D  microstructure. 
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Figure  5.7:  (a)  Macroscopic  stress-strain  response  by  plane  strain  and  plane  stress  VCFEM  simu¬ 
lation,  (b)  Stress-strain  hardening  rate  plots  for  the  Considere  condition. 
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Figure  5.8:  (a)  Experimental  micrographs,  (b)  VCFEM  simulated  micrograph  showing  damage  and 
contour  plot  of  effective  plastic  strain  at  8.88%  strain  in  section  1,  (c)  histogram  of  number  fraction 
of  cracked  particles  as  a  function  of  particle  size  by  Weibull  based  probabilistic  criterion,  and  (d) 
contour  plot  of  particle  fracture  probability  of  section  1  at  8.88%  strain. 
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Figure  5.9:  Number  fraction  of  cracked  particles  as  a  function  of  straining. 


5.6  Characteristic  Size  of  Microstructures 


The  influence  region  of  local  morphology  on  the  mechanical  response  is  characterized  by  a  mi- 
crostructural  representative  material  element  (RME)  that  is  critical  in  delineating  length  scales. 
The  RME  depicts  a  region  which  is  assumed  to  be  representative  of  the  entire  microstructure. 
Functions  that  distinguish  between  variations  in  stress/strain  distributions  for  local  disturbances 
in  microstructural  patterns  can  provide  important  insight  on  microstructure-property  relations. 
Marked  correlation  functions,  discussed  in  [22,  10,  35,  26]  for  multivariate  characterization  of  pat¬ 
terns,  are  evaluated  to  characterize  length  scales  or  RME  size  in  the  presence  of  damage.  A  mark 
may  be  identifled  with  an  appropriate  microstructural  variables,  e.g.  in  this  case  a  variable  that 
related  to  quantification  of  damage.  The  marked  correlation  function  for  a  heterogeneous  domain 
W  of  volume  V  containing  N  heterogeneities  is  mathematically  expressed  as  [22,  23,  10]: 


M(r)  = 


s(t)  ’ 


1  V 


N  k' 


mimk{r)  and  g{r) 


»=i  fc=i 


1  dKjr) 
47rr2  dr 


(5.8) 


where  K{r)  is  the  second  order  intensity  function  defined  in  [34,  35],  g{r)  is  the  pair  distribution 
function  and  H{r)  is  the  mark  intensity  function.  The  H{r)  function  reduces  to  the  K{r)  function  if 
all  heterogeneities  have  the  same  mark.  A  mark  associated  with  the  heterogeneity  is  denoted  as 
rui,  k'  is  the  number  of  heterogeneities  which  have  their  centers  within  a  sphere  of  radius  r  around 
heterogeneity,  for  which  the  mark  is  mk,  and  m  is  the  mean  of  all  marks.  By  definition  M{r) 
establishes  a  relation  between  the  location  and  associated  variables  for  heterogeneities.  Two  marks 
are  considered  in  this  study.  The  first  corresponds  to  particle  cracks  and  are  designated  as  m,-  =  1 
for  a  cracked  particle  and  mi  =  2  for  an  intact  particle.  The  second  corresponds  to  the  probability 
of  particle  fracture,  which  signifies  the  propensity  to  advance  the  microstructural  damage  state. 


99 


Figure  5.10:  Simulations  showing  effective  plastic  strain(%)  and  cracked  particles  in  three  different 
subsets  of  the  entire  micrograph  of  section  1,  (a)  RME  0  with  dimension  195fim  x  155/im  (b)  RME 
1  with  dimension  150/xm  x  155/im  and  (c)  RME  2  with  dimension  116/xrna:115^m. 


The  Af(r)  function  statistically  stabilizes  at  near-unit  values  at  a  distance  r,„ter  at  which  the  local 
morphology  ceases  to  have  any  significant  influence  on  evolving  variable.  Values  of  M(r)  >  1  show 
positive  correlation,  while  M(r)  <  1  indicates  repulsion  between  marks.  This  distance  r inter  is  an 
indication  of  the  physical  range  of  interaction  and  is  significant  in  making  decisions  about  length 
scales  and  the  RME  size. 

The  marked  correlation  functions  corresponding  to  cracked  particles  and  probability  of  fracture, 
are  plotted  in  figures  ??a  and  b  from  the  simulation  of  the  entire  micrograph  of  section  1  (termed 
as  RME  0)  with  dimensions  195/im  x  ISS/um.  The  dotted  line  corresponds  to  the  unit  M{t) 
for  uniform  distribution  of  spherical  heterogeneities  with  identical  marks.  Contour  plots  of  the 
equivalent  plastic  strain  in  the  simulated  micrograph  with  cracked  particles  are  shown  in  figure 
5.11a.  The  particle  area  fraction  for  this  micrograph  18.37%  and  the  total  number  of  particles 
and  cracked  particles  are  105  and  34  respectively.  The  plots  are  made  with  only  upto  40%  of 
the  entire  micrograph,  or  80/j.m  to  avoid  boundary  effects  in  M{r).  The  M{r)  functions  in  both 
figures  approximately  stabilize  at  near-unit  values  at  a  distance  r,nter  of  about  60  fim.  At  this 
distance,  the  local  morphology  is  expected  to  have  a  significantly  reduced  influence  on  the  evolving 
variables.  The  slower  attenuation  of  M(r)  for  particle  fracture  at  shorter  range  indicates  the  strong 
effect  of  the  local  morphology  on  damage  evolution.  Next,  a  smaller  region  (RME  1)  is  selected 
for  damage  simulation  corresponding  to  the  stabilized  region  in  the  M{r)  plots.  Since  the  stable 
region  is  60^m,  the  dimension  of  the  micrograph  is  chosen  to  be  150/im  x  Ibbfim,  incorporating 
the  scaling  factor,  i.e.  ^  =  150.  This  is  shown  with  the  box  in  figure  5.11a.  Again  the  contour 
plots  of  plastic  strain  with  cracked  particles  by  VCFEM  simulation  are  shown  in  figure  5.11b.  The 
dominant  crack  behavior  is  quite  similar  to  that  for  RME  0,  even  though  there  is  some  difference 
near  the  boundary.  Also  the  plastic  strain  contours  and  limiting  values  are  similar.  The  particle 
area  fraction  for  RME  1  is  18.13%  with  the  total  number  of  particles  and  cracked  particles  at  84 
and  28  respectively.  The  M(r)  plots  in  figure  ??c  and  d  show  that  the  functions  may  still  be 
assumed  to  stabilize  at  around  60/im.  A  smaller  subset  (RME  2),  with  dimensions  116/ima:115/im 
is  next  simulated  and  the  plastic  strain  is  depicted  in  figure  5.11c.  Significantly  different  plastic 
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strains  and  cracking  pattern  is  observed  for  this  microstructure.  The  particle  area  fraction  for 
this  micrograph  is  18.68%  with  a  total  of  51  particles  of  which  18  are  cracked.  The  plots  of  M(r) 
function  do  not  stabilize  in  the  domain  of  the  simulation  window.  Through  this  analysis  the  size 
effect  of  microstructure,  needed  for  adequate  representation  and  analysis  in  the  presence  of  evolving 
damage  is  demonstrated. 

5.7  Conclusions 

In  this  work,  a  combination  of  experimental  and  computational  methods  are  utilized  to  characterize 
and  understand  the  evolution  of  microscopic  damage  that  cause  failure  in  naturally  aged  commercial 
SiC  particle  reinforced  DRA’s.  The  main  mode  of  damage  for  the  naturally  aged  material  is  found 
to  be  particle  cracking.  Larger  particles  in  particle  rich  regions  are  more  susceptible  to  cracking 
than  those  in  particle  sparse  regions.  Spatial  distribution  of  particles  plays  a  more  important  role 
in  damage  than  particle  size  for  this  material.  A  sensitivity  analysis  with  respect  to  microstructural 
parameters  infers  that  the  strongest  influence  on  particle  cracking  comes  from  the  size  and  local 
spatial  distribution.  Particle  shape,  orientation  and  nearest  neighbor  orientation  have  relatively 
smaller  effect  on  damage  initiation.  Histograms  of  particles  forming  the  dominant  damage  path  in 
comparison  with  aU  cracked  particles  reveal  that  larger  particles  oriented  in  loading  direction  and  in 
relatively  rich  areas  are  more  susceptible  to  contribute  to  a  dominant  crack  in  the  microstructure. 
In  an  attempt  to  identify  discriminating  characteristics  of  2D  micrographs  that  may  be  of  helpful 
in  making  dominant  damage  predictions  for  the  actual  3D  microstructures,  probability  density 
functions  of  particle  size  nearest  neighbor  distance,  and  second  order  intensity  function  K{r)  of 
spatial  distribution  are  plotted.  Better  representation  of  damage  is  possible  with  those  sections 
that  have  higher  peaks  at  lower  near  neighbor  distances  and  longer  tails,  as  well  as  have  propensity 
towards  clustering. 

Next  the  two  dimensional  Voronoi  cell  finite  element  model  is  used  to  simulate  microstructural 
damage  evolution  in  computer  generated  equivalent  micrographs.  Both  macroscopic  and  micro¬ 
scopic  variables  obtained  by  the  VCFEM  simulation  are  compared  with  experimental  observations. 
The  macroscopic  stress-strain  plot  for  the  plane  strain  analysis  is  found  to  yield  quite  good  match 
with  experiments  if  the  difference  in  the  initial  yield  strength  due  to  plastic  constraint  is  subtracted. 
Prediction  of  the  onset  of  plastic  instability  by  the  Considere  criterion  is  also  found  to  be  in  reason¬ 
ably  good  agreement  with  the  experimental  results.  For  the  microstructural  results  with  number 
of  cracked  particles  in  different  size  ranges,  the  Weibull  model  is  found  to  give  better  concurrence 
with  experiments.  A  plot  of  the  number  fraction  of  cracked  particles  as  a  function  of  strairung 
shows  that  at  lower  strains  sections  with  particle  rich  regions  damage  rapidly,  but  the  rate  slows 
down  with  additional  deformation.  Finally,  the  marked  correlation  functions  are  evaluated  to  char¬ 
acterize  length  scales  and  representative  material  element  size  in  the  presence  of  damage.  Particle 
cracks  and  the  probability  of  particle  fracture  are  chosen  to  be  the  marks.  The  study  reveals  that  a 
significantly  large  portion  of  the  microstructure  should  be  analyzed  for  reasonable  accuracy  in  the 
presence  of  damage.  The  correlation  functions  do  not  stabilize  below  a  certain  length  scales  and 
this  keeps  growing  with  increased  damage.  In  summary,  various  important  characteristics  of  grow¬ 
ing  damage  are  investigated  in  this  work  to  understand  the  role  of  microstructure  in  the  material 
failure  process. 
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Figure  5.11:  Marked  correlation  function  as  a  function  of  radial  distance;  (a)  cracked  particles  as 
marks  and  (b)  probability  of  cracking  as  makr  for  RME  0;  (c)  cracked  particles  as  marks  and  (d) 
probability  of  cracking  as  makr  for  RME  1;  (e)  cracked  particles  as  marks  and  (f)  probability  of 
cracking  as  makr  for  RME  2. 
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